Questions on QSVT algorithm and Pennylane implementation

[jag]

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.

[CatalinaAlbornoz]

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.

[Guillermo_Alonso]

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

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

[jag]

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.

[Guillermo_Alonso]

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

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

[jag]

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]]

[CatalinaAlbornoz]

Hi @jag ,

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

[CatalinaAlbornoz]

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.

[jag]

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.

[CatalinaAlbornoz]

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!

[jag]

Hi @CatalinaAlbornoz ,

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

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

@CatalinaAlbornoz

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.

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]?

---

@CatalinaAlbornoz

If you do this you get the right eigenvalues in the imaginary part.

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

@CatalinaAlbornoz

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.

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.

[Danial_Motlagh]

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

[jag]

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.

[CatalinaAlbornoz]

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

[jag]

@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]]

[CatalinaAlbornoz]

Hi @jag ,

I’m not sure if this will work but how about using something like the sum_even_odd_circ function in the QSVT in practice demo?

And thanks for sharing the working code for the odd degree polynomial.

[jag]

Hi @CatalinaAlbornoz ,

Welcome! I hope the code would be useful for somebody.

I tried the approach you stated from the QSVT in Practice | PennyLane Demos demo for even function func = lambda x: 1/2 * (1.0 + 0.0*x + 1.0*np.square(x)) but the output is not as expected. The approach in the demo seems to be for better approximation. I’m not super sure. Maybe I’m doing something wrongly.

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

np.set_printoptions(threshold=np.inf)

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=', '))

poly = PolyTaylorSeries().taylor_series(
    func=func,
    degree=11,
    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, "proj_anc":1, "anc": 1})

block_encoding = qml.BlockEncode(A, wires=reg["b"].tolist())
phi1, phi2 = phi_qsvt[: len(phi_qsvt) // 2], phi_qsvt[len(phi_qsvt) // 2 :]

dim = A.shape[0] if A.ndim > 0 else 1
projectors_even = [qml.PCPhase(angle, dim= dim, wires=reg["b"].tolist()) for angle in phi1]
projectors_odd = [qml.PCPhase(angle, dim= dim, wires=reg["b"].tolist()) for angle in phi2]

def circ():
    qml.Hadamard(wires=reg["proj_anc"])  # equal superposition

    qml.ctrl(qml.QSVT, control=(reg["proj_anc"][0],), control_values=(0,))(block_encoding, projectors_even)
    qml.ctrl(qml.QSVT, control=(reg["proj_anc"][0],), control_values=(0,))(block_encoding, projectors_odd)

    qml.Hadamard(wires=reg["proj_anc"])  # equal superposition

    # These gates applied to the ancilla is to encode the 
    # polynomial transformation in the real part of the output matrix.
    # Alternatively, one can multiple the matrix to be block-encoded by i
    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["proj_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=', '))

[jag]

@CatalinaAlbornoz

For debugging, I removed all the code related to pyqsp generation of angles to narrow down the problem. I’m using the higher level Pennylane function qml.qsvt — PennyLane 0.42.0 documentation.

The classically computed output from SVD is equivalent to the qsvt output for odd degree function but doesn’t work for even degree function. I’m not sure whether this is a bug in the implementation or something wrong in my code.

import scipy
import scipy.linalg
import pennylane as qml
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]]

# DOESN'T WORK: Even degree function
coeff = [0.5, 0.0, 0.5]

# WORKS: Odd degree functions
# coeff = [0.0, 0.5, 0.0, 0.5]

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

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

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

def circ():
    qml.qsvt(A, coeff, encoding_wires=reg["b"].tolist(), block_encoding="embedding")
    return qml.state()

U_A = qml.matrix(circ, wire_order=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=', '))

[CatalinaAlbornoz]

Hi @jag ,

Thank you for adding this new example. I think the issue is in how you’re calculating the expected output because if we look at the last example in the qml.qsvt docs your expected output doesn’t match the results. I’m not exactly sure what’s wrong, I’m looking into it. I’ll let you know what I can find.

[CatalinaAlbornoz]

Hi @jag ,

I was able to chat with my colleague Danial and he explained why the comparison wasn’t working for even polynomials!

QSVT was designed to work with non-square matrices. If you want to get a polynomial of a matrix A with dimensions (m,n) you can’t really do A^2=A \times A because you can’t multiply (m,n) \times (m,n). Since this typical definition for A^2 doesn’t work, they define a new one within the QSVT framework. We know that A=WSV where W, S and V correspond to U, s and Vh when you use scipy.linalg.svd. Hence A^2 is defined (in the QSVT framework) as: A^2=(V^\dagger S W^\dagger)WSV. If you look at this closely you’ll notice that W^\dagger and W cancel out, and SS=S^2. So A^2=(V^\dagger S W^\dagger)WSV=V^\dagger S SV=V^\dagger S^2V. This is the key that makes QSVT work, but it’s also why the expected output for even and odd polynomials is different.

Here are some examples of how this works for different polynomials:

• A=WSV
• A^2=(V^\dagger S W^\dagger)WSV=V^\dagger S^2V
• A^3=WSV(V^\dagger S W^\dagger)WSV=(WSV)V^\dagger S^2V=WS^3V
• A^4=(V^\dagger S W^\dagger)WSV(V^\dagger S W^\dagger)WSV=(V^\dagger S W^\dagger)WS^3V=V^\dagger S^4V

These interleaving sets of (WSV) and (V^\dagger S W^\dagger) are the blocks that you will often see in examples and papers.

It also means that the expected output for an odd polynomial is W Poly(S) V, while for an even polynomial the output is V^\dagger Poly(S) V.

Hence you need to compare with a different expected output
This is an example of how you could test this with code:

# Define your polynomial
# P(x) = -0.1 + 0.2 x^2 + 0.5 x^4
poly_even = np.array([-0.1, 0, 0.2, 0, 0.5])
# P(x) = -0.1x + 0.2 x^3 + 0.5 x^5
poly_odd = np.array([0, -0.1, 0, 0.2, 0, 0.5])

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

# Define your matrix
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]])

# Create your qsvt circuit
dev = qml.device("default.qubit")

@qml.qnode(dev)
def circuit():
    qml.qsvt(A, poly, encoding_wires=[0, 1, 2, 3, 4], block_encoding="embedding")
    return qml.state()

matrix = qml.matrix(circuit, wire_order=[0, 1, 2, 3, 4])()

# Find the ground truth
W, s, Vh = scipy.linalg.svd(A)
Vh_dagger = np.conjugate(Vh.T)

# Make your comparisons
print("SVD singular values:", s)
if parity == "odd":
  print("Expected output Odd:\n", np.array2string(np.round(W @ scipy.linalg.diagsvd(np.poly1d(poly[::-1])(s), *A.shape) @ Vh,3), separator=', '))
else:
  print("Expected output Even:\n", np.array2string(np.round(Vh_dagger @ scipy.linalg.diagsvd(np.poly1d(poly[::-1])(s), *A.shape) @ Vh,3), separator=', '))

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

If you’re working with Hermitian (hence square) matrices you can actually use the traditional definition for A^2 by using GQSP, which is very powerful and state of the art . We have a GQSP function in PennyLane which you can use. My colleague Danial in fact is the author on the paper on GQSP. The review section in that paper has a lot of math that might be useful to understand that method. If you have any questions on this please let us know.

I hope this helps!