Questions on QSVT algorithm and Pennylane implementation

Hi,

I have a couple questions on QSVT algorithm and Pennylane’s implementation. Referring to Theorem 26 from the paper Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics by Gilyén et al.

What I understood is that QSVT creates a matrix which is a close approximation of the singular value decomposition, i.e \tilde \Pi U \Pi \approx W {\cdot} P(\Sigma) {\cdot} V^\dagger where P(\Sigma) is the singular values after the polynomial function is applied to them. So, I created a short code to experiment as below. The polynomial function that I’m applying is P(x)=x.

import numpy as np
import pennylane as qml
import scipy

gen = np.random.default_rng(42)
A = gen.uniform(-1, 1, (4, 4))  # random 4x4 matrix
# target poly: P(x) = x
target_poly = [0,1]
wire_order = list(range(5))

U, S, Vh = scipy.linalg.svd(A)

classical_output = U @ np.diag([np.polyval(target_poly[::-1], s) for s in S]) @ Vh
print("Classical output:\n", classical_output)

norm = scipy.linalg.norm(A, ord=2)

if norm > 1:
    A_normalized = A / norm

print("\nnorm of A:", norm)

quantum_output = qml.matrix(qml.qsvt, wire_order=wire_order)(
    A_normalized, target_poly, encoding_wires=wire_order, block_encoding="embedding"
) # block-encoded in 5-qubit system

print("\nQuantum_output:\n", quantum_output[:4, :4].real*norm)

However, the output differs between the classical and quantum computation.

  1. Is this the expected output? My expectation is that the quantum output will be a close approximation of W {\cdot} P(\Sigma) {\cdot} V^\dagger. I’m not sure whether my understanding is correct or whether my code is wrong.

  2. When I tried with a polynomial with an imaginary coefficient, e.g. target_poly = [0,1j], I get an error AssertionError: Array must not have an imaginary part. Is this by design of the implementation in Pennylane or is this because QSVT algorithm only works with real coefficients?

Looking forward to your kind reply. Thank you.

Hi @jag,

According to this paper the QSVT algorithm only works for polynomials with real coefficients and definite parity. So this limitation would be general and not specific to PennyLane’s implementation.

Regarding the first question let me get back to you.

Hi @jag ! Taking a look to the documentation I saw that QSVT works if A is hermitian
Here is an updated script to check that it is actually working :grinning_face_with_smiling_eyes:

import numpy as np
import pennylane as qml
import scipy

gen = np.random.default_rng(42)
A = gen.normal(size=(4, 4)) 
A = (A + A.T) / 2 # generating hermitian matrix

target_poly = [0,1]
wire_order = list(range(5))

U, S, Vh = scipy.linalg.svd(A)

classical_output = U @ np.diag([np.polyval(target_poly[::-1], s) for s in S]) @ Vh
print("Classical output:\n", classical_output)

norm = scipy.linalg.norm(A, ord=1)

if norm > 1:
    A_normalized = A / norm

print("\nnorm of A:", norm)

quantum_output = qml.matrix(qml.qsvt, wire_order=wire_order)(
    A_normalized, target_poly, encoding_wires=wire_order, block_encoding="embedding"
) # block-encoded in 5-qubit system

print("\nQuantum_output:\n", quantum_output[:4, :4].real*norm)
Classical output:
 [[ 0.30471708 -1.49550965  0.36682502  0.50329771]
 [-1.49550965 -1.30217951 -0.36260176  0.40549931]
 [ 0.36682502 -0.36260176  0.87939797  0.62265064]
 [ 0.50329771  0.40549931  0.62265064 -0.85929246]]

norm of A: 3.5657902238248567

Quantum_output:
 [[ 0.30471708 -1.49550965  0.36682502  0.50329771]
 [-1.49550965 -1.30217951 -0.36260176  0.40549931]
 [ 0.36682502 -0.36260176  0.87939797  0.62265064]
 [ 0.50329771  0.40549931  0.62265064 -0.85929246]]

Also note that you have to use the norm-1 instead of the norm-2

2 Likes

Hi @CatalinaAlbornoz and @Guillermo_Alonso,

Thank you for the helpful paper reference and the code. The code works now. I have a couple of follow-up questions which I encountered while experimenting with the QSVT algorithm and it’s implementation in Pennylane to understand it in practice.

  1. When I tried a constant polynomial, e.g. P(x)=1 where I replace the target_poly with target_poly = [1], I get an error, i.e. AssertionError: The polynomial must have at least degree 1. In similar veins as my previous questions, is this a algorithmic limitation of QSVT or a design choice of the implementation in Pennylane?
  2. For complex matrixes, the line to generate Hermitian matrix A = (A + A.T) / 2 will be replaced with A = (A + A.conj().T) / 2, right?

Thank you.

The function P(x) = 1 does not make use of block encoding—it’s equivalent to applying no operator at all. So in that case, applying QSVT doesn’t make sense :grinning_face_with_smiling_eyes:

On the other hand, yes, for complex matrices it would be that, but it’ll be harder to interpret the resulting matrix

1 Like

Hi @Guillermo_Alonso, thank you for the response. If I understood correctly, QSVT transforms the singular values of a matrix, i.e. \tilde \Pi U \Pi \approx W {\cdot} P(\Sigma) {\cdot} V^\dagger. However, I couldn’t figure out why the code below differs between the resultant matrix after singular value transformation using classical SVD computation and using QSVT? I’m trying to debug it for long but couldn’t find the issue. Is it because pyqsp in this case is using Chebyshev basis or other reasons?

import scipy
import scipy.linalg
import pennylane as qml
from pyqsp import angle_sequence
from pyqsp.poly import PolyTaylorSeries
import numpy as np

A = [[ 0.15421977, -0.66132613,  0.307225  , -0.66669095],
 [-0.98743463, -0.07939858,  0.0695699 , -0.11759568],
 [-0.0085032 , -0.4568372 ,  0.54552881,  0.70258507],
 [ 0.03341867,  0.58961292,  0.77664093, -0.21924566]]
func = lambda x: 1.0 + 0.0*x + 1.0*np.square(x)

A = A / np.linalg.norm(A, ord=1)
W, s, Vh = scipy.linalg.svd(A)

print("SVD singular values:", s)
print("Expected output:\n", np.array2string(np.round(W @ scipy.linalg.diagsvd(func(s), *A.shape) @ Vh,3), separator=', '))

poly = PolyTaylorSeries().taylor_series(
    func=func,
    degree=100,
    max_scale=1.0,
    chebyshev_basis=True,
    cheb_samples=200)

(phi_pyqsp, red_phiset, parity) = angle_sequence.QuantumSignalProcessingPhases(
    poly,
    signal_operator="Wx",
    method='sym_qsp',
    chebyshev_basis=True,
    )

phi_qsvt = qml.transform_angles(phi_pyqsp, "QSP", "QSVT")  # convert pyqsp angles to be compatible with QSVT

wire_order = list(range(3))
block_encoding = qml.BlockEncode(A, wires=wire_order)
print("Block encoded matrix:\n", np.array2string(block_encoding.matrix()[:4, :4], separator=', '))
projectors = [qml.PCPhase(angle, dim=A.shape[0], wires=wire_order) for angle in phi_qsvt]
U_A = qml.matrix(qml.QSVT, wire_order=wire_order)(block_encoding, projectors)

print("QSVT transformed matrix (real part):\n", np.array2string(np.round(U_A[:4, :4].real, 3), separator=', '))
print("QSVT transformed matrix (imaginary part):\n", np.array2string(np.round(U_A[:4, :4].imag, 3), separator=', '))

Output:

Expected output:
 [[ 0.203, -0.868,  0.403, -0.875],
 [-1.297, -0.104,  0.091, -0.154],
 [-0.011, -0.6  ,  0.716,  0.923],
 [ 0.044,  0.774,  1.02 , -0.288]]

Block encoded matrix:
 [[ 0.08629249, -0.37003998,  0.1719054 , -0.37304182],
 [-0.55251149, -0.04442687,  0.03892731, -0.06579976],
 [-0.0047579 , -0.25561976,  0.30524647,  0.3931261 ],
 [ 0.01869916,  0.3299134 ,  0.43456349, -0.12267723]]

QSVT transformed matrix (real part):
 [[-0.508, -0.   , -0.   , -0.   ],
 [-0.   , -0.508, -0.   , -0.   ],
 [-0.   , -0.   , -0.508,  0.   ],
 [-0.   , -0.   ,  0.   , -0.508]]

QSVT transformed matrix (imaginary part):
 [[ 0.657,  0.   ,  0.   ,  0.   ],
 [ 0.   ,  0.657,  0.   ,  0.   ],
 [ 0.   ,  0.   ,  0.657, -0.   ],
 [ 0.   ,  0.   , -0.   ,  0.657]]

Hi @jag ,

I’m not sure but someone in our team will get back to you on this next week.

Hi @jag ,

Sorry for the delay in our response. It seems that there may be an issue in how you’re calculating the coefficients, however it’s hard to tell. One of my colleagues will be taking a deeper look next week to see if we can find anything specific. We’ll let you know if we can find it.

Hi @CatalinaAlbornoz ,

Thank you for getting back to me. I’m trying hard to figure out what is the issue also. Looking forward to your assistance.

Have a pleasant weekend. :sunflower:

Hi @jag ,

My colleague Danial figured out the issue.

On one hand, the polynomial should be bounded with norm less or equal than 1 in the interval [-1,1]. So you need to multiply the polynomial (func) by 1/2.
If you do this you get the right eigenvalues in the imaginary part.

Even degree polynomials don’t give the same “shape” but the singular values are the same. So in these cases it might be best to use qml.BlockEncode() and then qml.GQSP() so that the matrix matches the way you would expect.

Let me know if this works for you!

Hi @CatalinaAlbornoz ,

Thank you for getting back to me. I’m really lost here. :grimacing:

I start off with sharing my code and then below it, I highlight what I’ve done and my questions.

import scipy
import scipy.linalg
import pennylane as qml
from pyqsp import angle_sequence
from pyqsp.poly import PolyTaylorSeries
import numpy as np

A = [[ 0.15421977, -0.66132613,  0.307225  , -0.66669095],
 [-0.98743463, -0.07939858,  0.0695699 , -0.11759568],
 [-0.0085032 , -0.4568372 ,  0.54552881,  0.70258507],
 [ 0.03341867,  0.58961292,  0.77664093, -0.21924566]]
func = lambda x: 1/2 * (1.0 + 0.0*x + 1.0*np.square(x))

A = A / np.linalg.norm(A, ord=1)
W, s, Vh = scipy.linalg.svd(A)

print("SVD singular values:", s)
print("Expected output:\n", np.array2string(np.round(W @ scipy.linalg.diagsvd(func(s), *A.shape) @ Vh,3), separator=', '))
print("Eigenvalues of matrix:\n", np.array2string(scipy.linalg.eigvals(A), separator=","))

poly = PolyTaylorSeries().taylor_series(
    func=func,
    degree=10,
    max_scale=1.0,
    chebyshev_basis=True,
    cheb_samples=200)

(phi_pyqsp, red_phiset, parity) = angle_sequence.QuantumSignalProcessingPhases(
    poly,
    signal_operator="Wx",
    method='sym_qsp',
    chebyshev_basis=True,
    )

phi_qsvt = qml.transform_angles(phi_pyqsp, "QSP", "QSVT")  # convert pyqsp angles to be compatible with QSVT

wire_order = list(range(3))
block_encoding = qml.BlockEncode(A, wires=wire_order)
print("Block encoded matrix:\n", np.array2string(block_encoding.matrix()[:4, :4], separator=', '))
projectors = [qml.PCPhase(angle, dim=A.shape[0], wires=wire_order) for angle in phi_qsvt]
U_A = qml.matrix(qml.QSVT, wire_order=wire_order)(block_encoding, projectors)

print("QSVT transformed matrix (real part):\n", np.array2string(np.round(U_A[:4, :4].real, 3), separator=', '))
print("QSVT transformed matrix (imaginary part):\n", np.array2string(np.round(U_A[:4, :4].imag, 3), separator=', '))

First part

I multiplied func by 1/2, which resulted in func = lambda x: 1/2 * (1.0 + 0.0*x + 1.0*np.square(x)).

Question: Isn’t Chebyshev basis which is used in the pyqsp part of the code is already bounded by [-1,1]?


In my code, I’ve computed the eigenvalues of the matrix classically

print("Eigenvalues of matrix:\n", scipy.linalg.eigvals(A))

which gives the eigenvalues

Eigenvalues of matrix:
 [-0.44594767+0.33796201j,-0.44594767-0.33796201j, 0.5581651 +0.03923467j,
  0.5581651 -0.03923467j]

Question: Which doesn’t seem to match the imaginary part of the QSVT computed matrix as below?

QSVT transformed matrix (imaginary part):
 [[ 0.657,  0.   ,  0.   ,  0.   ],
 [ 0.   ,  0.657,  0.   ,  0.   ],
 [ 0.   ,  0.   ,  0.657, -0.   ],
 [ 0.   ,  0.   , -0.   ,  0.657]]

Second part

Question: From what I understand, I keep using the qml.BlockEncode() as per my code?


I read the documentation of qml.GQSP() and it says

  • 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.

Question: How can I reshape my pyqsp generated angles to be compatible with qml.GQSP?


Question: What if I have an odd degree polynomial?


Goal

Just to reiterate, what I’m trying to see is from the code, the Expected output matrix that I computed classically is (approximately) equal to the real part of the matrix computed by QSVT.

Why are you comparing with the eigenvalues of the original matrix? You should be comparing with the singular values of the original matrix transformed by your polynomial.

If you want the Expected output to match the matrix computed just choose an odd polynomial like 0.5x + 0.5x^3 and look at the imaginary part of the matrix outputted by QSVT.

Whether the function is encoded in the real or imaginary part really depends on the block-encoding convention being used. You can always switch between them by multiplying the transformed block by i which would correspond to applying X S X on the ancilla qubit after the QSP sequence

Thank you for the response. I got confused on comparing with eigenvalues as the previous answer mentioned about eigenvalue. Nevertheless, I got it working for odd polynomials, the Expected output is same as computed by QSVT.

Question: Is it possible for the Expected output to be the same for even polynomials, e.g. 1/2 * (1.0 + 0.0*x + 1.0*np.square(x)) using qml.QSVT? Catalina mentioned about using qml.GQSP() but the angle required by GQSP seems to be of shape (3, d+1) which is not how pyqsp library generates angles. Hence, I’m unsure how to use qml.GQSP.

I’m doing these experiments because I’m cross-checking whether the polynomial that I’m using is correctly applied to the matrix with QSVT especially at values close to 0. The polynomials that I’ve stated in my code are toy examples.

Thanks for adding the additional details @jag !

I’ll check to see if there’s something we can do about even polynomials. I’m glad that at least you got it working for odd polynomials. If you want to add your solution here it might help others in the future (and it might help us look into alternatives for even polynomials).

@CatalinaAlbornoz Thank you for your continued support.

For anyone encountering the same questions posed by me in this thread, below is the code that works for odd degree polynomial where the classically computed matrix with transformed singular values is (approximately) equal to the real part of the complex matrix produced by QSVT.

import scipy
import scipy.linalg
import pennylane as qml
from pyqsp import angle_sequence
from pyqsp.poly import PolyTaylorSeries
import numpy as np

A = [[ 0.15421977, -0.66132613,  0.307225  , -0.66669095],
 [-0.98743463, -0.07939858,  0.0695699 , -0.11759568],
 [-0.0085032 , -0.4568372 ,  0.54552881,  0.70258507],
 [ 0.03341867,  0.58961292,  0.77664093, -0.21924566]]

# Odd degree polynomials
func = lambda x: 0.5*x + 0.5*(x**3)
# func = lambda x: x

A = A / np.linalg.norm(A, ord=1)
W, s, Vh = scipy.linalg.svd(A)

print("SVD singular values:", s)
print("Expected output:\n", np.array2string(np.round(W @ scipy.linalg.diagsvd(func(s), *A.shape) @ Vh,3), separator=', '))

poly = PolyTaylorSeries().taylor_series(
    func=func,
    degree=10,
    max_scale=1.0,
    chebyshev_basis=True,
    cheb_samples=200)

(phi_pyqsp, red_phiset, parity) = angle_sequence.QuantumSignalProcessingPhases(
    poly,
    signal_operator="Wx",
    method='sym_qsp',
    chebyshev_basis=True,
    )

phi_qsvt = qml.transform_angles(phi_pyqsp, "QSP", "QSVT")  # convert pyqsp angles to be compatible with QSVT

reg = qml.registers({"b": 3, "anc": 1})

wire_order = list(range(3))
block_encoding = qml.BlockEncode(A, wires=reg["b"].tolist())
projectors = [qml.PCPhase(angle, dim=A.shape[0], wires=reg["b"].tolist()) for angle in phi_qsvt]

def circ():
    qml.QSVT(block_encoding, projectors)

    # These gates applied to the ancilla is to encode the 
    # polynomial transformation in the real part of the output matrix.
    # Alternatively, one can multiply the matrix to be block-encoded by i, i.e. 1j in Python
    qml.PauliX(wires=reg["anc"])
    qml.adjoint(qml.S)(wires=reg["anc"])
    qml.PauliX(wires=reg["anc"])

    return qml.state()

U_A = qml.matrix(circ, wire_order=reg["anc"].tolist() + reg["b"].tolist())()

print("QSVT transformed matrix (real part):\n", np.array2string(np.round(U_A[:4, :4].real, 3), separator=', '))
print("QSVT transformed matrix (imaginary part):\n", np.array2string(np.round(U_A[:4, :4].imag, 3), separator=', '))

Output:

SVD singular values: [0.55954235 0.55954235 0.55954235 0.55954234]
Expected output:
 [[ 0.057, -0.243,  0.113, -0.245],
 [-0.363, -0.029,  0.026, -0.043],
 [-0.003, -0.168,  0.2  ,  0.258],
 [ 0.012,  0.217,  0.285, -0.081]]
QSVT transformed matrix (real part):
 [[ 0.057, -0.243,  0.113, -0.245],
 [-0.363, -0.029,  0.026, -0.043],
 [-0.003, -0.168,  0.2  ,  0.258],
 [ 0.012,  0.217,  0.285, -0.081]]
QSVT transformed matrix (imaginary part):
 [[ 0.141, -0.603,  0.28 , -0.607],
 [-0.9  , -0.072,  0.063, -0.107],
 [-0.008, -0.416,  0.497,  0.64 ],
 [ 0.03 ,  0.537,  0.708, -0.2  ]]

However, even degree polynomials, e.g. func = lambda x: 1/2 * (1.0 + 0.0*x + 1.0*np.square(x)) doesn’t produce the expected output using QSVT. @CatalinaAlbornoz, it will be helpful if someone have some idea how to make it work for even degree polynomials.

Output:

SVD singular values: [0.55954235 0.55954235 0.55954235 0.55954234]
Expected output:
 [[ 0.101, -0.434,  0.202, -0.438],
 [-0.648, -0.052,  0.046, -0.077],
 [-0.006, -0.3  ,  0.358,  0.461],
 [ 0.022,  0.387,  0.51 , -0.144]]
QSVT transformed matrix (real part):
 [[ 0.657,  0.   ,  0.   ,  0.   ],
 [ 0.   ,  0.657,  0.   ,  0.   ],
 [ 0.   ,  0.   ,  0.657, -0.   ],
 [ 0.   ,  0.   , -0.   ,  0.657]]
QSVT transformed matrix (imaginary part):
 [[ 0.486,  0.   ,  0.   ,  0.   ],
 [ 0.   ,  0.486,  0.   ,  0.   ],
 [ 0.   ,  0.   ,  0.486, -0.   ],
 [ 0.   ,  0.   , -0.   ,  0.486]]