Quantum Chebyshev Transform

What are Chebyshev polynomials?

The \(n\) -th order Chebyshev polynomial of the first kind is defined as

Chebyshev polynomials of the first kind are a set of orthogonal polynomials. They are complete on the interval \([-1,1]\), meaning that any polynomial \(f(x)\) of degree \(N\) can be written exactly as a series in \(T_n(x)\) up to order \(N\) as \(f(x) = \sum_{j=0}^N a_j T_j(x)\) on that interval.

We can write out the first few orders explicitly.

The recursion relation in the last line can be used to compute the next orders. Observe that odd and even order \(T_n\) are odd and even functions, respectively. The roots of the \(T_n(x)\) occur at the values

These are known as the Chebyshev nodes.

The nodes are plotted above along with the corresponding polynomials. Note that the polynomials are normalized such that $T_n(1)=1$, and they satisfy the following discrete orthogonality condition on the nodes of \(T_N(x)\) for \(k, \ell<N\)

The Chebyshev polynomials have a lot of nice properties that make them a good choice for expanding a function into a series on a real interval. They are often utilized in numerical integration and interpolation methods, where the discrete Chebyshev transformation is needed to map between the function evaluated on a grid of points and the coefficients of the Chebyshev expansion of that function. Let’s take a look at the discrete transform.

Discrete Chebyshev transform

A common definition of the discrete Chebyshev transform uses the grid of Chebyshev nodes, also known as the roots grid, as the sampling points. The transform computes the coefficients for the \((N-1)\) -degree expansion of a function \(f(x)\) sampled at the \(N\) nodes of \(T_N(x)\). That is, the \(j\) -th coefficient, for \(0<j<N\), is

and for \(j=0\), there is a slightly different normalization factor

The inverse of the transform is then just the series expansion evaluated on the grid

Since a function expanded in Chebyshev polynomials in this way will be sampled on the non-uniformly spaced Chebyshev nodes, it will have more resolution at the boundary than in the middle. In fact, the grid of Chebyshev nodes minimizes the problem of Runge’s phenomenon. This can be beneficial if you are, for example, solving a differential equation and expect more interesting features at the boundary. Furthermore, the Chebyshev expansion gives the best polynomial approximation to a continuous function under the maximum norm.

In general, computing the expansion of a function in a complete set on a classical computer for a discrete number of sampling points would take \(\mathcal{O}(N^2)\) operations 🐌. However, because of the way the Chebyshev polynomials are defined in terms of cosine, the discrete Chebyshev transformation is related to the discrete cosine transform. This allows the discrete Chebyshev transform to be implemented in a way that leverages the efficiency of the fast-Fourier-transform-style algorithms for expansion, which take \(\mathcal{O}(N \log N)\) operations 🚀.

We can see the relation to the cosine transform by plugging in the definition of the Chebyshev polynomials and the nodes into the inverse transform and simplifying. Starting with the polynomials

then, using the definition of the nodes we obtain

Finally, we can use the cyclical property of cosine to convert a \(j \pi\) term in the argument to a \((-1)^{j}\) factor in the coefficient, resulting in

which looks just like a discrete cosine transform.

The quantum analogue of the discrete Chebyshev transform, the quantum Chebyshev transform, inherits the relation to the Fourier transform, allowing the transform to be designed efficiently by utilizing the quantum Fourier transform. Next we will discuss the quantum Chebyshev basis, where the Chebyshev polynomials appear in the state amplitudes.

Quantum Chebyshev basis

We can define the \(j\) -th Chebyshev basis state using \(N\) qubits as

where \(|k\rangle\) are the computational basis states and \(x_j^\mathrm{Ch}\) is the \(j\) -th node of the Chebyshev polynomial of order \(2^N-1\). Notice how the amplitudes of the basis state components are the Chebyshev polynomials evaluated at the \(j\) -th Chebyshev node, and that the normalization of the \(|0\rangle\) component is different from the rest, like in the classical transform. Due to the orthogonality of the Chebyshev polynomials and the normalization factors used, this construction guarantees the states are orthonormal, that is

The quantum Chebyshev transform circuit described in Ref. 1 maps computational basis states \(\{|x_j\rangle\}_{j=0}^{2^N-1}\) to Chebyshev basis states \(\{|\tau(x_j^\mathrm{Ch})\rangle\}_{j=0}^{2^N-1}\). Our goal is to design a circuit in PennyLane that applies the operation \(\mathcal{U}_\mathrm{QChT} = \sum_{j=0}^{2^N-1} |\tau(x_j^\mathrm{Ch})\rangle\langle x_j|\).

Designing the transform circuit

Let’s start from the end and look at the circuit diagram generated from the code we want to write. An auxiliary qubit is required, which will be the \(0\) indexed qubit, and the rest compose the state \(|x\rangle\) which starts in the computational basis, shown below as \(|\psi\rangle\). We demonstrate for \(N=4\) non-auxiliary qubits.

The intuition for the structure of the above circuit comes from the link between the discrete Chebyshev transform and the discrete cosine transform. Notice the use of the quantum Fourier transform (QFT) applied on all qubits. The quantum Chebyshev transform is an extended QFT circuit with some added interference and mixing of the elements. Note the auxiliary qubit starts and ends in the state \(|0\rangle\), and the amplitudes of the transformed state are all real valued. Let’s break down the circuit above into pieces that we will use inside our circuit function.

First, a Hadamard gate is applied to the auxiliary qubit, and then a CNOT ladder is applied, controlled on the auxiliary qubit. To start, we will define a function for the CNOT ladder.

import pennylane as qml

# Number of qubits (non-auxiliary qubit).
N = 4


def CNOT_ladder():
    for wire in range(1, N + 1):
        qml.CNOT([0, wire])

After the initial CNOT ladder comes an \(N+1\) QFT circuit, which can be implemented using qml.QFT.

Next are phase rotations and shifts. For the auxiliary qubit, there is a \(Z\) rotation by and angle of \(-\pi(2^N - 1)/2^{N+1}\) followed by a phase shift of \(-\pi/2^{(N+1)}\) . The other qubits are rotated in \(Z\) by \(\pi/2^{(j+1)}\), where \(j\) is the index of the qubit as labelled in the circuit diagram.

import numpy as np

pi = np.pi


def rotate_phases():
    """Rotates and shifts the phase of the auxiliary qubit and rotates the jth qubit by
    pi/2^(j+1) in Z."""
    qml.RZ(-pi * (2**N - 1) / 2 ** (N + 1), wires=0)
    qml.PhaseShift(-pi / 2 ** (N + 1), wires=0)
    for wire in range(1, N + 1):
        qml.RZ(pi / 2 ** (wire + 1), wires=wire)

Now a permutation of the qubits is used to reorder them. This is built using a multicontrolled NOT gate applied to each qubit from the initial state, which is controlled on the auxiliary qubit and all qubits with larger index than the target. The multicontrolled NOT gate can be implemented using a multicontrolled Pauli X gate. Let’s see what that looks like.

def permute_elements():
    """Reorders amplitudes of the conditioned states."""
    for wire in reversed(range(1, N + 1)):
        control_wires = [0] + list(range(wire + 1, N + 1))
        qml.MultiControlledX(wires=(*control_wires, wire))

In the above code, we use reversed to loop over the qubits in reverse order, to apply the controlled gate to the last qubit first. After the permutation is another CNOT ladder, which we already have a function for.

The last part is a phase adjustment of the auxiliary qubit: a rotation in \(Y\) by \(\pi/2\), a phase shift of \(-\pi/2\) and a multicontrolled \(X\) rotation by \(\pi/2\). All of the other qubits control the \(X\) rotation, but the control is sandwiched by Pauli \(X\) operators. We can implement the multicontrolled \(X\) rotation by using the function qml.ctrl on qml.RX, specifying the target wire in qml.RX and the control wires as the second argument of qml.ctrl.

def adjust_phases():
    """Adjusts the phase of the auxiliary qubit."""
    qml.RY(-pi / 2, wires=0)
    qml.PhaseShift(-pi / 2, wires=0)

    # First Pauli X gates.
    for wire in range(1, N + 1):
        qml.PauliX(wires=wire)

    # Controlled RX gate.
    qml.ctrl(qml.RX(pi / 2, wires=0), range(1, N + 1))

    # Second Pauli X gates.
    for wire in range(1, N + 1):
        qml.PauliX(wires=wire)

All together, we can construct the circuit. We have added qml.BasisState to initialize the input in any computational basis state with the optional argument state.

def QChT():
    """Performs the quantum Chebyshev transform."""
    qml.Hadamard(wires=0)
    CNOT_ladder()
    qml.QFT(wires=range(N + 1))
    rotate_phases()
    permute_elements()
    CNOT_ladder()
    adjust_phases()


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


@qml.qnode(dev)
def circuit(state=0):
    qml.BasisState(state=state, wires=range(1, N + 1))
    QChT()
    return qml.state()

Finally, we can reproduce the circuit diagram shown at the beginning of this section using qml.draw_mpl.

def circuit_to_draw():
    qml.BasisState(state=0, wires=range(1, N + 1))
    QChT()


fig, ax = qml.draw_mpl(circuit_to_draw, decimals=2, style="pennylane")()
fig.show()

Note we defined a new function for the circuit to simplify the drawing, removing the returned qml.state.

Testing the quantum Chebyshev transform

With our quantum Chebyshev transform circuit, let’s first check if the auxiliary qubit ends in the state \(|0\rangle\) and the output state amplitudes are real valued. To do this, we’ll input the computational basis state \(|7\rangle\), which will transform into \(|\tau(x_7^\mathrm{Ch})\rangle\). We expect the full output state to be \(|0\rangle|\tau(x_7^\mathrm{Ch})\rangle\), which means the second half of the amplitude vector should be zero (since there should be no state components with the auxiliary qubit in \(|1\rangle\)).

j = 7  # Initial state in computational basis.

total_state = circuit(state=j)  # State with auxiliary qubit.

# Round very small values to zero.
total_state = np.where(np.abs(total_state) < 1e-12, 0, total_state)
print(np.round(total_state, 3))

[ 0.25 -0.j  0.035-0.j -0.347+0.j -0.103+0.j  0.327-0.j  0.167-0.j
 -0.294-0.j -0.224+0.j  0.25 +0.j  0.273+0.j -0.196-0.j -0.312+0.j
  0.135+0.j  0.338+0.j -0.069-0.j -0.352+0.j  0.   +0.j  0.   +0.j
  0.   +0.j  0.   +0.j  0.   +0.j  0.   +0.j  0.   +0.j  0.   +0.j
  0.   +0.j  0.   +0.j  0.   +0.j  0.   +0.j  0.   +0.j  0.   +0.j
  0.   +0.j  0.   +0.j]

Indeed, we see the second half of the amplitude vector is zero. Furthermore, the first \(2^N\) entries are real valued, but let’s check if the amplitudes of the state components in the computational basis agree with our definition. Remember we expect

We can just check the circuit output by computing each state component directly.

def T_n(x, n):
    """Chebyshev polynomial of order n."""
    return np.cos(n * np.arccos(x))


def ch_node(j):
    """The jth node of Chebyshev polynomial 2^N - 1."""
    return np.cos(pi * (2 * j + 1) / 2 ** (N + 1))


def tau_amplitudes(x, k):
    """Computes the state amplitudes of tau in the computational basis."""
    prefactor = np.where(k == 0, 1 / 2 ** (N / 2), 1 / 2 ** ((N - 1) / 2))
    return prefactor * T_n(x, k)


# Reduce state size, effectively removing the auxiliary qubit.
state = np.real(total_state[: 2**N])  # Discard the small imaginary components.

k = range(2**N)  # Computational basis indices.
expected_state = tau_amplitudes(ch_node(j), np.array(k))

assert np.allclose(expected_state, state, atol=1e-9)  # Compare circuit output to calculated values.

The output state from the circuit is exactly what we want.

To give some intuition for how the quantum Chebyshev transform is related to the quantum Fourier transform, we can plot the overlap of \(|\tau(x_7^\mathrm{Ch})\rangle\) with the computational basis states \(|k\rangle\). This should look like a cosine plot, since Chebyshev polynomials are just the cosine of an argument which is proportional to \(k\).

import matplotlib.pyplot as plt

plt.style.use("pennylane.drawer.plot")

ks = np.linspace(0, max(k), 100)  # Continuous list of points.

fig = plt.figure(figsize=(6.4, 3.2))
ax = fig.add_axes((0.15, 0.23, 0.80, 0.72))  # Make room for caption.
ax.plot(k, state, "o", label="circuit")
ax.plot(ks, [tau_amplitudes(ch_node(j), kk) for kk in ks], label="interpolation")
ax.set(xlabel=r"$|k\rangle$", ylabel=r"Overlap $\langle k|\tau(x_7^\mathrm{Ch})\rangle$")
ax.legend()
fig.text(0.5, 0.05,
    r"Figure 4. Overlaps of $|\tau(x_7^\mathrm{Ch})\rangle$ with computational basis states.",
    horizontalalignment="center",
    size="small", weight="normal",
)
plt.show()

Notice that, unlike the quantum Fourier transform, these amplitudes are real instead of complex values. Also note that the overlap at \(|0\rangle\) is discontinuous because the \(|0\rangle\) amplitude of the Chebyshev state was adjusted in the definition to guarantee orthonormality of the Chebyshev basis.

Next, let’s test that orthonormality by computing the overlap at the nodes with all other Chebyshev basis states \(|\tau(x_j^\mathrm{Ch})\rangle\).

# Compute the overlap with the other basis states using np.vdot().
js = list(range(int(len(state))))
overlaps = [np.vdot(state, circuit(state=i)[: 2**N]) for i in js]

We can compare these circuit-calculated overlaps to the definition, plotting the squared overlaps at all values of \(x\). This continuous overlap function can be derived analytically as 1

where \(\tau(x)\) is a generalization of one of Chebyshev basis states defined earlier and \(x\) can be any value in \([-1,1]\) rather than just one of the nodes.

def overlap_sq(x, xp):
    """Computes the squared overlap between Chebyshev states."""
    numerator = T_n(xp, 2**N + 1) * T_n(x, 2**N) - T_n(xp, 2**N) * T_n(x, 2**N + 1)
    return numerator**2 / (2 ** (2 * N)) / (xp - x) ** 2

fig = plt.figure(figsize=(6.4, 3.2))
ax = fig.add_axes((0.15, 0.25, 0.8, 0.6))  # Make room for the caption.
ax.set(xlabel=r"x", ylabel="Squared Overlap")
ax_top = ax.twiny()

# Plot squared overlaps computed in the circuit.
nodes = [ch_node(i) for i in js]
ax.plot(nodes, np.abs(overlaps) ** 2, marker="o", label="circuit")

# Plot expected squared overlaps.
xs = np.linspace(-1, 1, 1000)
ax.plot(xs, [overlap_sq(x, nodes[j]) for x in xs], label="expectation")

# Set the second axis to show computational basis state indices.
ax_top.set_xlim(ax.get_xlim())
tick_labels = [str(i) for i in js]
tick_labels[0::2] = [''] * len(tick_labels[0::2])  # Omit even-numbered labels.
ax_top.set(xlabel=r"Chebyshev node $j$",
           xticks=nodes, xticklabels=tick_labels)
ax_top.grid(False)

ax.legend()
fig.text(0.5, 0.05,
    "Figure 5. Squared overlap of Chebyshev states "
    + r"$|\langle \tau(x_k^\mathrm{Ch})|\tau(x_7^\mathrm{Ch})\rangle|^2$",
    horizontalalignment="center",
    size="small", weight="normal",
)
plt.show()

We can see that the squared overlap between the basis states and the \(j=7\) state \(|\tau(x_7^\mathrm{Ch})\rangle\) is 0, unless \(x=x_7^\mathrm{Ch}\approx 0.1\), then the overlap is 1.

Conclusion

In this tutorial, we’ve gone through how to implement the quantum Chebyshev transform from the paper by Williams et al., and tested the circuit output by looking at the state amplitudes and the orthonormality. The properties of Chebyshev polynomials and the speed at which the quantum Chebyshev transform can be implemented make the Chebyshev basis an interesting choice of state space for quantum algorithms, such as generative modelling of probability distributions. To build a generative model in the Chebyshev basis, one could implement the quantum Chebyshev feature map from the same paper 1, which prepares a state in the Chevyshev space via a parameter \(x\).

References

1(1,2,3,4)

Chelsea A. Williams, Annie E. Paine, Hsin-Yu Wu, Vincent E. Elfving and Oleksandr Kyriienk. “Quantum Chebyshev transform: mapping, embedding, learning and sampling distributions.” arxiv:2306.17026 (2023).

2

Jorge J. Martínez de Lejarza, Hsin-Yu Wu, Oleksandr Kyriienko, Germán Rodrigo, Michele Grossi. “Quantum Chebyshev probabilistic models for fragmentation functions.” arxiv:2503.16073 (2025).