Questions on QSVT algorithm and Pennylane implementation

Hi @CatalinaAlbornoz ,

Thank you very much! I was really lost. :sweat_smile: It is way clearer now.

I have a couple of questions. My Pennylane version is 0.42.0-dev52.

Question 1: I’m using qml.GQSP for an arbitrary polynomial that I took from the documentation. However, the output that I’m computing using the classical SVD using scipy doesn’t tally with the output from GQSP? I’ve checked whether the matrix is Hermitian and it is Hermitian. My code is as below.

import numpy as np
import pennylane as qml
import scipy

# Qiskit is just used for visualization
from qiskit.visualization import array_to_latex

A = np.array([[-0.1, 0, 0, 0.1], [0, 0.2, 0, 0], [0, 0, -0.2, -0.2], [0.1, 0, -0.2, -0.1]])

print("Is A Hermitian:", scipy.linalg.ishermitian(A))

# P(x) = 0.1 + 0.2j x + 0.3 x^2
poly = [0.1, 0.2j, 0.3]

angles = qml.poly_to_angles(poly, "GQSP")

dev = qml.device("default.qubit")

@qml.prod # transforms the qfunc into an Operator
def unitary(wires):
    qml.BlockEncode(A, wires=wires)

@qml.qnode(dev)
def circuit():
    qml.GQSP(unitary(wires=[0,1,2]), angles, control = 3)
    return qml.state()

qml.draw_mpl(circuit)()

W, s, Vh = scipy.linalg.svd(A)

expected_output = np.round(W @ scipy.linalg.diagsvd(np.poly1d(poly[::-1])(s), *A.shape) @ Vh,3)
print("Expected output real:\n")
display(array_to_latex(expected_output.real, max_size=10000))
print("Expected output imaginary:\n")
display(array_to_latex(expected_output.imag, max_size=10000))

matrix = qml.matrix(circuit, wire_order=[0,1,2,3])()
print("GQSP output real\n")
display(array_to_latex(np.round(matrix[:4, :4].real, 3), max_size=10000))
print("GQSP output imaginary\n")
display(array_to_latex(np.round(matrix[:4, :4].imag, 3), max_size=10000))

Output:

Question 2: You mentioned that if I use Hermitian (square matrix), then I can use GQSP.

For an arbitrary matrix A (not necessarily Hermitian), can I do Hermitian dilation as following, block encode H using BlockEncode or Fabel and then use it with GQSP?

H = \begin{bmatrix} 0 & A \\ A^{\dagger } & 0 \end{bmatrix}

Question 3: Why doesn’t qml.poly_to_angles allow a degree 0 polynomial such as P(x) = 1x^0? If I run the following code:

import pennylane as qml

poly = [1]

angles = qml.poly_to_angles(poly, "GQSP")

there will be an error AssertionError: The polynomial must have at least degree 1. I’m thinking of a use case to do a polynomial transformation like WPoly(S)V = WIV = WV where I is the identity matrix and Poly(S)=1.

Question 4: As far as I know, pyqsp doesn’t have the ability to generate angles for GQSP which expects

  • angles (tensor[ float ]) – array of angles defining the polynomial transformation. The shape of the array must be (3, d+1), where d is the degree of the polynomial.

Also, pyqsp doesn’t accept function with complex coefficients such as func = lambda x: 0.1 + 0.2j*x + 0.3*x**2

On the other hand, Pennylane doesn’t have a feature to generate polynomial for an arbitrary function unless I generate phase angles from optimization as done in QSVT in Practice or if I already know the coefficients of the function.

Hence, just restricting to real coefficients, I am generating the angles using the following code:

import pennylane as qml
from pyqsp.poly import PolyTaylorSeries

# Some arbitrary function
# pyqsp doesn't accept complex coefficients
# func = lambda x: 0.1 + 0.2j*x + 0.3*x**2
func = lambda x: 0.1 + 0.2*x + 0.3*x**2

poly, scale = PolyTaylorSeries().taylor_series(
    func=func,
    degree=11, # EVEN degree like 10 doesn't work 
    max_scale=1.0,
    chebyshev_basis=True,
    ensure_bounded=True,
    return_scale=True,
    cheb_samples=200)

# Pennylane expects polynomial in standard basis - https://discuss.pennylane.ai/t/error-when-using-qml-qsvt/8612/5?u=jag
target_poly = np.polynomial.chebyshev.cheb2poly(poly.coef.tolist())

angles = qml.poly_to_angles(target_poly, "GQSP")

Is this an acceptable way to generate the angles to be used with GQSP?

Also, poly_to_angles throw an error ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (2,) + inhomogeneous part. when I use an even number, e.g. 10, as the degree argument value in the code above.

Thank you again for your kind help.

Hi @jag ,

  1. My colleague Guillermo mentioned that if you apply a polynomial over a block encoding with GQSP then you obtain the polynomial applied over the generic block encoding, which is not the same. This is a bit cryptic but I think the point is that BlockEncode has a very specific block encoding, and applying GQSP is basically doing a block encoding of the block encoding, which messes things up. I’ll try to get more clarity on this and get back to you.
  2. I think this might result in the same issues as question 1 so I’ll have to check on that.
  3. I think this may have been by design. I can see in the GQSP paper that “d = 0 is trivial” (eq 8) so some of the equations may just not be defined for degree 0.
  4. Can you post the full error you get for this one? I want to check where it’s originating.

I hope to come back with more answers soon!

Hi @CatalinaAlbornoz ,

Thank you for following up on my questions.

  1. & 2. I see in the documentation of GQSP that it expects a unitary as the input.

Parameters:

  • unitary (Operator) – the operator to be encoded by the GQSP circuit

Without block or fabel encoding used in conjunction with Hermitian dilation or some similar methods, how can an arbitrary (not necessarily Hermitian) matrix, if possible, can be used with GQSP? I’ll wait for more clarity from you.

  1. I read back the paper. If I understand correctly, to obtain the transformation WPoly(S)V = WIV = WV where I is the identity matrix and Poly(S) = 1, it will be an application of a rotation on the first qubit (Fig 2. from the paper below) without the need for further operations, i.e. I find \phi and \theta such that e^{i(\lambda + \phi)}cos(\theta)I = 1 in equation 8 of the paper? Hence, it’s a trivial operation that doesn’t need GQSP? The reason I’m conjecturing this is because I couldn’t conclude whether GQSP or QSVT doesn’t work for constant polynomial transformation, e.g. P(x)=1.

  1. The error is as attached. I’m using Pennylane version 0.42.0-dev52.

Hi @jag ,

I think we’re back to expecting the wrong expected output (plus a wire ordering detail).

What’s the true expected output?
We know that applying GQSP to a unitary U creates a matrix where the top left is Poly(U) as shown in the docs. So in this case we shouldn’t be expecting Poly(A) in the top left (as you had in your previous code), we’re expecting Poly(U).

Toy example:
We created U via BlockEncode, so we know that
\begin{align} U(A) &= \begin{bmatrix} A & \sqrt{I-AA^\dagger} \\ \sqrt{I-A^\dagger A} & -A^\dagger \end{bmatrix}. \end{align}

Since A is Hermitian (in this case), if we do the math for U^2 we’ll see that U^2=I.
So we can write a toy example just using a polynomial that is 0.99x^2

# P(x) = 0. + 0. x + 0.99 x^2

If we apply GQSP(U) then we expect the output to contain Poly(U) in the top left, which means we expect 0.99*U^2=0.99I.

On the other hand if our polynomial was # P(x) = 0. + 0.99 x + 0. x^2 then we would expect to get 0.99U on the top left.

Why this matters
This matters because, while GQSP works well for both odd and even polynomials, the way you compute your expected output differs. You again go back to the same situation as before where the expected output for an odd polynomial is WPoly(S)V, while for an even polynomial the output is V^\dagger Poly(S)V.

The code below should illustrate this:

import numpy as np
import pennylane as qml
import scipy

# Qiskit is just used for visualization
from qiskit.visualization import array_to_latex

# Define our Hermitian matrix an print it
A = np.array([[-0.1, 0, 0, 0.1], [0, 0.2, 0, 0], [0, 0, -0.2, -0.2], [0.1, 0, -0.2, -0.1]])

print("A\n",np.round(A,3))
print("Is A Hermitian:", scipy.linalg.ishermitian(A))

# BlockEncode A to create a unitary that can be used with GQPE
U = qml.BlockEncode(A, wires=[1,2,3])
print("Is U unitary: ", qml.is_unitary(U))
print("Is U Hermitian: ", qml.is_hermitian(U))

# Find the matrix for the unitary and print it
u = U.matrix(wire_order=[1,2,3])
print("u\n",np.round(u,3))

# Define your polynomial
# P(x) = 0. + 0. x + 0.99 x^2
poly_even = [0.0, 0., 0.99]
# P(x) = 0. + 0.99 x + 0. x^2
poly_odd = [0.0, 0.99, 0.]

# Select your polynomial
parity = "odd" # choose between "odd" or "even"
if parity == "odd":
  poly = poly_odd
  print("Odd polynomial")
else:
  poly = poly_even
  print("Even polynomial")

# Calculate the angles for GQSP
angles = qml.poly_to_angles(poly, "GQSP")

# Do GQSP on the unitary
dev = qml.device("default.qubit")

@qml.qnode(dev)
def circuit():
    qml.GQSP(U, angles, control = 0)
    return qml.state()

# Calculate the expected output
W, s, Vh = scipy.linalg.svd(u)
Vh_dagger = np.conjugate(Vh.T)
if parity == "odd":
  expected_output_u = np.round(W @ scipy.linalg.diagsvd(np.poly1d(poly[::-1])(s), *u.shape) @ Vh,3)
  print("Expected output Odd:\n")
  display(array_to_latex(expected_output_u.real, max_size=10000))
else:
  expected_output_u = np.round(Vh_dagger @ scipy.linalg.diagsvd(np.poly1d(poly[::-1])(s), *u.shape) @ Vh,3)
  print("Expected output Even:\n")
  display(array_to_latex(expected_output_u.real, max_size=10000))

matrix = qml.matrix(circuit, wire_order=[0,1,2,3])()
print("GQSP output real\n")
display(array_to_latex(np.round(matrix.real, 5), max_size=10000))

What’s the wire ordering issue?
I noticed that things look good when the control is in qubit zero but the ordering in the matrices gets mixed up when the control is in qubit 3. I’m sure there’s math to explain this but I didn’t get the chance to look into it so if it’s not critical the easiest is just to choose qubit 0 as the control in GQSP for now.

Going back to question 2
Given the math I mentioned above then you could just factor in the new transformation you need to do in case the original matrix is not Hermitian. Let me know if this makes sense, otherwise I can reword it.

Let me know if this solves your first two questions!

Hi @CatalinaAlbornoz ,

Yes, it’s clearer now on Question 1 & 2.

I did a small experimentation to see clearly if doing GQSP(U) has any correspondence to doing the same Poly(S) transformation on the original matrix A. The slightly updated code is as below:

### Code above is the same ###
# Calculate the expected output
W, s, Vh = scipy.linalg.svd(u)
Vh_dagger = np.conjugate(Vh.T)

WA, sA, VhA = scipy.linalg.svd(A)
VhA_dagger = np.conjugate(VhA.T)

if parity == "odd":
  expected_output_A = np.round(WA @ scipy.linalg.diagsvd(np.poly1d(poly[::-1])(sA), *A.shape) @ VhA,3)
  print("Expected output A Odd:\n")
  display(array_to_latex(expected_output_A.real, max_size=10000))
  expected_output_u = np.round(W @ scipy.linalg.diagsvd(np.poly1d(poly[::-1])(s), *u.shape) @ Vh,3)
  print("Expected output u Odd:\n")
  display(array_to_latex(expected_output_u.real, max_size=10000))
else:
  expected_output_A = np.round(VhA_dagger @ scipy.linalg.diagsvd(np.poly1d(poly[::-1])(sA), *A.shape) @ VhA,3)
  print("Expected output A Even:\n")
  display(array_to_latex(expected_output_A.real, max_size=10000))
  expected_output_u = np.round(Vh_dagger @ scipy.linalg.diagsvd(np.poly1d(poly[::-1])(s), *u.shape) @ Vh,3)
  print("Expected output u Even:\n")
  display(array_to_latex(expected_output_u.real, max_size=10000))

matrix = qml.matrix(circuit, wire_order=[0,1,2,3])()
print("GQSP output real\n")
display(array_to_latex(np.round(matrix.real, 5), max_size=10000))

It appears that for the odd polynomial, the top-left block after GQSP(U) is the same as transforming the singular values of A by Poly(S):

However, for the even polynomial, the top-left block after GQSP(U) is not the same as transforming the singular values of A by Poly(S):

It seems to me that for these toy odd and even polynomials, there is no advantage in using GQSP over QSVT as QSVT can also do even polynomial transformation and it works directly with A?

Also, does GQSP have the hard requirement that the original matrix is Hermitian which requires transformations like Hermitian dilation if the original matrix is non-Hermitian whereas in QSVT, I can directly apply it to a block-encoded non-Hermitian original matrix?

Looking forward to your thoughts on question 3 and 4 too. I’m learning a lot from these discussions. Thank you very much.

Hi @jag ,

It looks like for GQSP we should be comparing with the eigenvalues instead of the singular values (this is new info for me too). So the comparison code I shared before is probably not the best way. I’ll share some code in the next few days.

My understanding is that GQSP(U)=Poly(A)=WPoly(D)W^\dagger, where A=WDW^\dagger, and D is the matrix of eigenvalues.

The advantage of using GQSP as opposed to QSVT lies in the lifting of restrictions on the polynomial, such as the definite parity restriction. It does however impose specific restrictions on the matrix, so this is probably the consideration you want to have for using one or the other.

For GQSP you need to use a unitary U. If your matrix is already Hermitian then you can block encode it into U

\begin{align} U(A) &= \begin{bmatrix} A & i\sqrt{I-A^2} \\ i\sqrt{I-A^2} & -A \end{bmatrix}. \end{align}

If you don’t have a Hemitian matrix A then you can indeed use an encoding such as:

\begin{align} A &= \begin{bmatrix} 0 &B \\ B^\dagger & 0 \end{bmatrix}. \end{align}

Note that this will require double the number of qubits compared to having a matrix that’s already Hermitian.

Because of how everything’s constructed, applying Hadamards on either side of GQSP would give you some linear combination of the even and odd polynomials of B. So if for example our target polynomial is 0.99x^2 then we would have something like H (GQSP(U)) H\approx 1/2 \sum(Poly_{even}(B)+Poly_{odd}(B)).

I’ll try to get some code running to check all of this but I hope it answers some questions.

RE question 3: the implementation for qml.poly_to_angles when selecting GQSP relies on Q_poly = scale_factor * Polynomial.fromroots(inside_circle) (see here in the docs). This only works if there are roots, hence why there’s a check to ensure that the polynomial has at least degree 1.

Regarding the error you’re getting in question 4 it looks like it might be because of using cheb2poly instead of poly2cheb. I’m not too convinced that this is the issue but it’s something you can try (I haven’t got the chance to code this yet).

I hope this hasn’t been too confusing, I’m also discovering new things as we go.

Note: my colleague Danial was a bit confused about your original goals in using QSVT and GQSP. He mentioned that the interest is generally in having U=e^{iH} where A=cos(H). This has some nice properties when you define U in terms of A as above, since \begin{align} U^n &= \begin{bmatrix} e^{inH} &0 \\ 0 & e^{-inH} \end{bmatrix} & \approx \begin{bmatrix} cos(nH) & i \cdot sin(nH) \\ i \cdot sin(nH) & -cos(nH) \end{bmatrix} \end{align}

I don’t know if this can be useful information for you, or if you just care about the polynomial transformation.

Hi guys! :slight_smile: GQSP is a somewhat advanced algorithm, and it’s understandable that there’s confusion. I’m going to try to shed some light on it.

We need to see GQSP as an algorithm that transforms the unitary U, not A. If we apply a polynomial P with GQSP, we’re performing a block encoding of P(U) — we don’t get P(A). What usually happens is that another polynomial R ends up being applied to A — that is, we’re actually doing a block encoding of R(A).

So the question we need to ask is: can we choose the polynomial R by making a good choice of P?

Well, it turns out this is possible if we choose the block encoding carefully! If we use qml.BlockEncode, it won’t work because U² = Identity, since the information about A gets lost — it’s no longer encoded anywhere.

Here’s where the magic comes in: what happens if instead of qml.BlockEncode(wires = [1,2,3]), we use qml.BlockEncode(wires = [1,2,3]) @ qml.Z(1)? This would be equivalent to flipping the sign of the second column of the block encoding. I encourage you to make the change and, with pen and paper, start computing U, , , etc. You’ll see that in this case, the polynomial being applied to A is A, 2A² - I, 4A³ - 3A, etc. And this is nothing other than the Chebyshev polynomials!

Therefore, there is a relationship between P and R: the Chebyshev transform :rocket:

I’m leaving you a script so you can see how it works

import numpy as np
from numpy.polynomial.chebyshev import poly2cheb
import pennylane as qml
import scipy

def matrix_polynomial(A, poly):
    "aux function to calculate poly(A) that I'm using to check the solution"
    n = A.shape[0]
    result = poly[0] * np.eye(n)
    Apow = np.eye(n)
    for i in range(1, len(poly)):
        Apow = Apow @ A
        result += poly[i] * Apow
    return result

A = np.array([[-0.1, 0, 0, 0.1], [0, 0.2, 0, 0], [0, 0, -0.2, -0.2], [0.1, 0, -0.2, -0.1]])
U = qml.BlockEncode(A, wires=[1,2,3]) @ qml.Z(1)

print("block encoded A:")
block_encoded_A = qml.matrix(U)[:4,:4]
print(block_encoded_A)


poly = np.array([0.1, 0.4, 0.4])
cheb = poly2cheb(poly)

# Calculate the angles for GQSP
angles = qml.poly_to_angles(cheb, "GQSP")

# Do GQSP on the unitary
dev = qml.device("default.qubit")
@qml.qnode(dev)
def circuit():
    qml.GQSP(U, angles, control = 0)
    return qml.state()

print("target -> P(A):")
target = matrix_polynomial(block_encoded_A, poly).real
print(target)

matrix = qml.matrix(circuit, wire_order=[0,1,2,3])()
print("GQSP output real\n")
output = np.round(matrix.real, 5)[:4,:4]
print(output)

print(np.allclose(output, target))
block encoded A:
[[-0.1  0.   0.   0.1]
 [ 0.   0.2  0.   0. ]
 [ 0.   0.  -0.2 -0.2]
 [ 0.1  0.  -0.2 -0.1]]
target -> P(A):
[[ 0.068  0.    -0.008  0.032]
 [ 0.     0.196  0.     0.   ]
 [-0.008  0.     0.052 -0.056]
 [ 0.032  0.    -0.056  0.084]]
GQSP output real

[[ 0.068 -0.    -0.008  0.032]
 [-0.     0.196 -0.     0.   ]
 [-0.008 -0.     0.052 -0.056]
 [ 0.032  0.    -0.056  0.084]]
True

Note that I decided to use an arbitrary polynomial without parity, main advantage of GQSP :grinning_face_with_smiling_eyes:

1 Like

@CatalinaAlbornoz ,

No worries. I appreciate your kindness in taking time to look into my questions.

@Guillermo_Alonso

Thank you for taking time to help me. I experimented with your code. However, correct me if I’m wrong, in your code, there doesn’t seem to any singular value transformation but rather it does the application of matrix polynomial.

My goal is to apply a polynomial to the singular values of an arbitrary matrix (non-necessarily Hermitian) such that A = WSV transforms to H = WPoly(S)V. Question: Is GQSP still applicable for this goal or QSVT is the go-to algorithm for this case?

Question: Also, can GQSP or QSVT work with a constant polynomial function, e.g. P(x) = 1, i.e. transform all singular values to the constant, 1 in this case? I can’t find a conclusive evidence that these algorithms doesn’t work with constant polynomial function. Probably, they are trivial transformations that doesn’t need these algorithms at all but this is just conjecturing from my part. As I explained previously, when I run:

import pennylane as qml
poly = [1]
angles = qml.poly_to_angles(poly, "GQSP") # or qml.poly_to_angles(poly, "QSVT") 

I get an error AssertionError: The polynomial must have at least degree 1. Catalina rightly pointed:

So, to reiterate, is this a limitation of the poly_to_angles function in Pennylane or GQSP / QSVT doesn’t work with constant polynomial transformation?

Hi @jag, let me check and get back to you.

Hi @jag ,

The degree 1 issue is a limitation in PennyLane. Since a degree 0 polynomial just gives the identity then there’s nothing really to do or compute. GQSP/QSVT define this as the trivial case, so technically it is not a limitation of these methods. Since this polynomial is trivial it’s just not implemented in PennyLane.

Regarding your other question, if your goal is specifically to apply a polynomial to the singular values, and not to the matrix itself, then GQSP is not the tool for the task.
As I mentioned before this gives you a polynomial of the eigenvalues and not the singular values. I show this with the code below. So if you strictly want the polynomial of the singular values the it might be better to stick with QSVT, with the limitation that it entails in the Polynomial.

# Workflow
import numpy as np
from numpy.polynomial.chebyshev import poly2cheb
import pennylane as qml
import scipy

# -------- Workflow -----------

# Select your Hermitian matrix and polynomial
A = np.array([[-0.1, 0, 0, 0.1], [0, 0.2, 0, 0], [0, 0, -0.2, -0.2], [0.1, 0, -0.2, -0.1]])
poly = np.array([0.1, 0.4, 0.4])

# Block encode your matrix into a unitary (GQSP requires a Unitary)
U = qml.BlockEncode(A, wires=[1,2,3]) @ qml.Z(1) # Note the addition of Z here
block_encoded_A = qml.matrix(U)[:4,:4]
print("Block encoded A:\n", block_encoded_A)
# Check that Block encoded A is still Hermitian
print("\nIs block_encoded_A Hermitian?", scipy.linalg.ishermitian(block_encoded_A)) # Must be Hermitian

# Calculate the angles for GQSP
cheb = poly2cheb(poly)
angles = qml.poly_to_angles(cheb, "GQSP")

# Do GQSP on the unitary
dev = qml.device("default.qubit")
@qml.qnode(dev)
def circuit():
    qml.GQSP(U, angles, control = 0)
    return qml.state()

matrix = qml.matrix(circuit, wire_order=[0,1,2,3])()
output = np.round(matrix.real, 5)[:4,:4]
print("\nGQSP output real\n", output)

# -------- Checks -----------

# Get the eigendecomposition of block_encoded_A
eigenvalues, eigenvectors  = np.linalg.eig(block_encoded_A)
ev_dag = np.conjugate(eigenvectors.T)

# Optionally check that A=P@D@P_dag where P is the matrix of eigenvectors and D is the matrix of eigenvalues
# Dei = scipy.linalg.diagsvd(eigenvalues, *A.shape)
# print("P@D@P_dag: ", np.round(eigenvectors@Dei@ev_dag,3))

# Polynomial of the eigenvalues
Poly_ei = np.poly1d(poly[::-1])(eigenvalues)
D_poly_ei = scipy.linalg.diagsvd(Poly_ei, *A.shape)

# The expected output is Poly(A)=P@Poly(D)@P_dag 
expected_output = np.round(eigenvectors @ D_poly_ei @ ev_dag,3) # A=P@D@P_dag so Poly(A)=P@Poly(D)@P_dag
print("\nExpected output:\n", expected_output.real)

# Optionally check directly against P(A)
def matrix_polynomial(A, poly):
    "aux function to calculate poly(A) that I'm using to check the solution"
    n = A.shape[0]
    result = poly[0] * np.eye(n)
    Apow = np.eye(n)
    for i in range(1, len(poly)):
        Apow = Apow @ A
        result += poly[i] * Apow
    return result

print("\ntarget -> P(A):")
target = matrix_polynomial(block_encoded_A, poly).real
print(target)

print("\nAre the output and expected_output equal?",np.allclose(output, expected_output))

Let me know if you have any further questions. This has been a very interesting discussion.

Hi @CatalinaAlbornoz ,

I was away for a while. I greatly appreciate @Guillermo_Alonso and your answers. It’s way clearer for me now. I’ll experiment further and reach out should I have further questions.

1 Like