From block encodings to quantum signal processing and quantum singular value transformation#

Code at: https://github.com/qiboteam/qibo/tree/main/examples/qsp_qsvt.

This tutorial follows a matrix through the three Qibo building blocks that turn it into a polynomial transformation of its spectrum on a quantum computer:

Block encoding (block_encoding_circuit, and the standalone linear combination of unitaries (LCU), lcu_circuit) embeds a matrix \(A\) in a block of a unitary matrix \(U\).
Quantum signal processing (QSP, qsp_phases and qsp_circuit) applies a function to the eigenphases of a walk operator \(W\) built from \(U\), for instance to simulate the time evolution \(e^{-i t H}\) of a Hamiltonian \(H\).
Quantum singular value transformation (QSVT, qsvt_phases and qsvt_circuit) applies a polynomial directly to the singular values of the matrix that \(U\) encodes.

                                   A is Hermitian: qubitization
                             ┌─────────────────────────────────────┐
                             │  W = i R U           (Section 2)    │
                             │  qsp_phases ──► qsp_circuit         │
matrix A ──► block encoding  │  f(θ) on the eigenphases of W,      │
             U, α, a         │  e.g. e^{-i t H}                    │
(Section 1)                  └─────────────────────────────────────┘
             ┌─────────────────────────────────────┐
             │  any A          (Section 3)         │
             │  qsvt_phases ──► qsvt_circuit       │
             │  P(ς) on the singular values of A/α │
             └─────────────────────────────────────┘

Conventions. A circuit is a block encoding of \(A\) with normalization \(\alpha\) if, when its \(a\) auxiliary qubits start and end in \(|0\rangle\), it acts on the \(n\) system qubits as \(A / \alpha\). In terms of the unitary matrix \(U\) of the circuit,

\[\frac{A}{\alpha} = \left( \langle 0|^{\otimes a} \otimes \mathbb{1} \right) U \left( |0\rangle^{\otimes a} \otimes \mathbb{1} \right),\]

where \(\mathbb{1}\) is the identity on the system qubits. The auxiliary qubits are always the first qubits of the circuit, so the encoded matrix is the top-left \(2^{n} \times 2^{n}\) block of \(U\). Since a block of a unitary matrix cannot have a singular value larger than one, the normalization is at least the largest singular value of \(A\). For a broader overview of how QSP leads to QSVT, see Ref. [5].

Setup#

We use the numpy backend. The phase computation of QSP and QSVT needs linear algebra routines, such as polynomial root finding, that some backends (for instance qibojit) do not implement.

[ ]:

from scipy.special import jv

from qibo import Circuit, gates, get_backend, set_backend
from qibo.models import (
    block_encoding_circuit,
    lcu_circuit,
    qsp_circuit,
    qsp_phases,
    qsvt_circuit,
    qsvt_phases,
)

set_backend("numpy")
backend = get_backend()

1. Block encoding#

block_encoding_circuit(matrix, method=...) takes a square matrix \(A\), padded with zeros to dimension \(2^{n}\) if needed, and returns the circuit, the normalization \(\alpha\) and the number \(a\) of auxiliary qubits. There are two methods:

method="dense": with \(B = A / \alpha\), whose singular values are at most one, the circuit is a single gate that holds the unitary dilation

\[\begin{split}U = \begin{pmatrix} B & \sqrt{\mathbb{1} - B B^{\dagger}} \\ \sqrt{\mathbb{1} - B^{\dagger} B} & -B^{\dagger} \end{pmatrix},\end{split}\]

where \(B^{\dagger}\) is the conjugate transpose of \(B\) and the square roots are matrix square roots. It uses \(a = 1\) auxiliary qubit. If \(\alpha\) is not given, it is the largest singular value of \(A\) when this exceeds one, and one otherwise. - method="lcu": the matrix is decomposed into Pauli strings \(P_{i}\), i.e. tensor products of \(\mathbb{1}\), \(X\), \(Y\) and \(Z\), as \(A = \sum_{i} c_{i} P_{i}\), with complex coefficients \(c_{i}\). The circuit is \(\mathrm{PREPARE}_{L}^{\dagger} \, \mathrm{SELECT} \, \mathrm{PREPARE}_{R}\): the two \(\mathrm{PREPARE}\) circuits load the weights \(\sqrt{|c_{i}| / \alpha}\) (the right one also with the phases of \(c_{i}\)) on \(a = \lceil \log_{2} L \rceil\) auxiliary qubits, where \(L\) is the number of non-zero terms, and \(\mathrm{SELECT}\) applies \(P_{i}\) controlled on the auxiliary qubits being in the basis state \(|i\rangle\). The normalization is fixed to \(\alpha = \sum_{i} |c_{i}|\) [4].

As an example, we encode the Hamiltonian of two qubits

\[H = -0.4 \, Z \otimes Z + 0.3 \, X \otimes \mathbb{1} + 0.3 \, \mathbb{1} \otimes X,\]

and check that the top-left block of the circuit is \(H / \alpha\).

[ ]:

pauli_x, pauli_z = backend.matrices.X, backend.matrices.Z
identity = backend.matrices.I()

hamiltonian = (
    -0.4 * backend.kron(pauli_z, pauli_z)
    + 0.3 * backend.kron(pauli_x, identity)
    + 0.3 * backend.kron(identity, pauli_x)
)

for method in ("dense", "lcu"):
    circuit, alpha, nauxiliary = block_encoding_circuit(hamiltonian, method=method)
    block = circuit.unitary(backend)[:4, :4]
    print(
        f"{method}: {circuit.nqubits} qubits, {nauxiliary} auxiliary, "
        f"{len(circuit.queue)} gates, alpha = {alpha:.2f}, "
        f"block equals H / alpha: {backend.allclose(block, hamiltonian / alpha)}"
    )

dense: 3 qubits, 1 auxiliary, 1 gates, alpha = 1.00, block equals H / alpha: True
lcu: 4 qubits, 2 auxiliary, 28 gates, alpha = 1.00, block equals H / alpha: True

Both methods build the encoding from the full matrix, so their classical cost grows exponentially with the number of qubits. When the Pauli decomposition is already known, which is the usual case for Hamiltonians, the LCU circuit can be built directly with lcu_circuit, which takes a list of circuits \(U_{i}\) on the same qubits and their coefficients \(c_{i}\), and block-encodes \(\sum_{i} c_{i} U_{i} / \alpha\). Gates that are already controlled, such as a CNOT, are allowed in \(U_{i}\), but every \(U_{i}\) must be exact, including its global phase.

[ ]:

paulis = []
for operators in ([gates.Z(0), gates.Z(1)], [gates.X(0)], [gates.X(1)]):
    pauli = Circuit(2)
    pauli.add(operators)
    paulis.append(pauli)

circuit, alpha, nauxiliary = lcu_circuit(paulis, [-0.4, 0.3, 0.3])
block = circuit.unitary(backend)[:4, :4]
print(backend.allclose(block, hamiltonian / alpha), alpha, nauxiliary)

True 1.0 2

2. From the block encoding to QSP#

2.1 Qubitization: the walk operator#

QSP does not act on the block encoding \(U\) itself, but on a walk operator \(W\) whose eigenphases carry the spectrum of \(H\). Let \(\Pi = |0\rangle\langle 0|^{\otimes a} \otimes \mathbb{1}\) be the projector on the auxiliary qubits being in \(|0\rangle\), and \(R = 2 \Pi - \mathbb{1}\) the reflection about it. If \(H\) is Hermitian and \(U\) is a Hermitian block encoding of \(H / \alpha\), then

\[W = i \, R \, U\]

acts, on the two-dimensional space spanned by \(|0\rangle|\lambda\rangle\) and \(U|0\rangle|\lambda\rangle\) for each eigenvector \(|\lambda\rangle\) of \(H\) with eigenvalue \(\lambda\), as a rotation with eigenvalues \(e^{\pm i \theta}\), where

\[\sin\theta = \frac{\lambda}{\alpha}.\]

This construction is called qubitization [2]. Both methods of block_encoding_circuit return a Hermitian \(U\) when \(A\) is Hermitian. Qibo does not provide the walk as a function yet, so we build it here. Since \(R = -\left( X^{\otimes a} \, C^{a-1}Z \, X^{\otimes a} \right)\), with \(C^{a-1}Z\) a \(Z\) gate controlled on the other \(a - 1\) auxiliary qubits, the walk is \(U\), followed by that multi-controlled gate and a global phase \(-i\). The phase has to be a gate, because QSP applies \(W\) under control, which turns a global phase into a relative one.

[ ]:

def qubitized_walk(block_encoding: Circuit, nauxiliary: int) -> Circuit:
    """Walk operator W = i R U of a Hermitian block encoding U."""
    nqubits = block_encoding.nqubits
    auxiliary = list(range(nauxiliary))

    walk = Circuit(nqubits)
    walk.add(block_encoding.on_qubits(*range(nqubits)))
    # sign flip of the auxiliary state |0...0>
    walk.add(gates.X(qubit) for qubit in auxiliary)
    walk.add(gates.Z(auxiliary[-1]).controlled_by(*auxiliary[:-1]))
    walk.add(gates.X(qubit) for qubit in auxiliary)
    # global phase -i, as a gate
    walk.add(gates.Unitary(-1j * identity, 0))

    return walk

We check the claim for the dense block encoding of \(H\): \(U\) is Hermitian, and the eigenphases \(\theta\) of the walk satisfy \(\alpha \sin\theta = \lambda\) for the eigenvalues \(\lambda\) of \(H\).

[ ]:

circuit, alpha, nauxiliary = block_encoding_circuit(hamiltonian, method="dense")
unitary = circuit.unitary(backend)
walk = qubitized_walk(circuit, nauxiliary)

eigenphases = backend.angle(backend.eigvals(walk.unitary(backend)))
sines = sorted({round(float(value), 4) for value in alpha * backend.sin(eigenphases)})
eigenvalues = [round(float(value), 4) for value in backend.eigvalsh(hamiltonian)]

print("U is Hermitian:", backend.allclose(unitary, backend.dagger(unitary)))
print("alpha * sin(theta):", sines)
print("eigenvalues of H:  ", eigenvalues)

U is Hermitian: True
alpha * sin(theta): [-0.7211, -0.4, 0.4, 0.7211]
eigenvalues of H:   [-0.7211, -0.4, 0.4, 0.7211]

2.2 QSP phases for the time evolution#

Given real trigonometric series \(A(\theta) = \sum_{k=0}^{d} a_{k} \cos(k\theta)\) and \(C(\theta) = \sum_{k=1}^{d} c_{k} \sin(k\theta)\) with \(A(0) = 1\) and \(|A + iC| \le 1\), qsp_phases(cosine, sine) returns \(N = 2d\) phases for which qsp_circuit applies \(A(\theta) + i C(\theta)\) to every eigenphase \(\theta\) of the walk [1]. The circuit queries the walk \(N\) times, under control of one extra qubit, the QSP qubit.

To obtain \(e^{-i \tau \sin\theta}\), whose value on the eigenphases of the walk is \(e^{-i \tau \lambda / \alpha}\), we use the Jacobi-Anger expansion, with \(J_{k}\) the Bessel function of the first kind of order \(k\):

\[\cos(\tau \sin\theta) = J_{0}(\tau) + 2 \sum_{k \ge 1} J_{2k}(\tau) \cos(2k\theta), \qquad \sin(\tau \sin\theta) = 2 \sum_{k \ge 0} J_{2k+1}(\tau) \sin\left((2k+1)\theta\right).\]

The target is \(A + iC = \cos(\tau \sin\theta) - i \sin(\tau \sin\theta)\), so the sine coefficients are the negative of those of \(\sin(\tau \sin\theta)\). The Bessel functions decay super-exponentially for orders larger than \(\tau\), so the series is truncated at the smallest degree \(d\) whose tail is below the target error \(\epsilon\). Like the walk, this step is not yet part of Qibo’s public API.

[ ]:

def evolution_phases(tau: float, epsilon: float):
    """QSP phases of exp(-i tau sin(theta)) with error below epsilon."""
    # the tail beyond the next 64 orders is negligible
    degree = 1
    while (
        2 * sum(abs(jv(order, tau)) for order in range(degree + 1, degree + 65))
        > epsilon
    ):
        degree += 1

    orders = backend.arange(degree + 1)
    bessel = backend.cast(jv(backend.to_numpy(orders), tau), dtype="float64")
    cosine = backend.where(backend.mod(orders, 2) == 0, 2 * bessel, 0.0)
    cosine[0] = bessel[0]
    sine = backend.where(backend.mod(orders, 2) == 1, -2 * bessel, 0.0)[1:]

    return qsp_phases(cosine, sine)

2.3 The QSP circuit#

The block encoding is of \(H / \alpha\), so the QSP circuit with parameter \(\tau\) applies \(e^{-i \tau H / \alpha}\). To simulate the time \(t\) we therefore choose \(\tau = \alpha t\). The qubits of the circuit are the QSP qubit, then the qubits of the walk, i.e. the auxiliary qubits of the block encoding and the system. The evolution is the block in which the QSP qubit and the auxiliary qubits start and end in \(|0\rangle\), which we extract by reshaping the unitary matrix of the circuit.

[ ]:

time = 1.5
target = backend.expm(-1j * time * hamiltonian)

phases = evolution_phases(alpha * time, 1e-8)
qsp = qsp_circuit(walk, phases)

nsystem = circuit.nqubits - nauxiliary
shape = (2, 2**nauxiliary, 2**nsystem) * 2
block = backend.reshape(qsp.unitary(backend), shape)[0, 0, :, 0, 0, :]
error = backend.matrix_norm(block - target, order=2)

print(
    f"dense: {len(phases)} walk queries, {qsp.nqubits} qubits, "
    f"error below 1e-7: {error < 1e-7}"
)

dense: 20 walk queries, 4 qubits, error below 1e-7: True

2.4 QSP with the LCU block encoding#

The same pipeline works with method="lcu", which does not need the dense matrix of the walk. There is one caveat: qsp_circuit adds the QSP qubit as a control to every gate of the walk, and Qibo cannot add a control to a gate that is already controlled, such as the CNOT, Toffoli or multi-controlled gates of the LCU circuit. Until qsp_circuit handles them, the walk can be wrapped in a single gate, which gives up the gate-level structure of the LCU circuit but verifies the pipeline.

[ ]:

circuit, alpha, nauxiliary = block_encoding_circuit(hamiltonian, method="lcu")
unitary = circuit.unitary(backend)
walk = qubitized_walk(circuit, nauxiliary)

dense_walk = Circuit(walk.nqubits)
dense_walk.add(gates.Unitary(walk.unitary(backend), *range(walk.nqubits)))

phases = evolution_phases(alpha * time, 1e-8)
qsp = qsp_circuit(dense_walk, phases)

nsystem = circuit.nqubits - nauxiliary
shape = (2, 2**nauxiliary, 2**nsystem) * 2
block = backend.reshape(qsp.unitary(backend), shape)[0, 0, :, 0, 0, :]
error = backend.matrix_norm(block - target, order=2)

print("U is Hermitian:", backend.allclose(unitary, backend.dagger(unitary)))
print(
    f"lcu: {len(phases)} walk queries, {qsp.nqubits} qubits, "
    f"error below 1e-7: {error < 1e-7}"
)

U is Hermitian: True
lcu: 20 walk queries, 5 qubits, error below 1e-7: True

3. From the block encoding to QSVT#

QSVT does not need a Hermitian matrix or a walk operator. Let \(A / \alpha = \sum_{i} \varsigma_{i} |w_{i}\rangle\langle v_{i}|\) be the singular value decomposition (SVD) of the encoded matrix, with singular values \(\varsigma_{i}\) and left and right singular vectors \(|w_{i}\rangle\) and \(|v_{i}\rangle\). For a real polynomial \(P\) of degree \(d\) with the parity of \(d\) and \(|P(x)| \le 1\) for \(x \in [-1, 1]\), qsvt_circuit applies

\[\begin{split}P^{(SV)}\left(\frac{A}{\alpha}\right) = \begin{cases} \sum_{i} P(\varsigma_{i}) \, |w_{i}\rangle\langle v_{i}| & d \text{ odd,} \\ \sum_{i} P(\varsigma_{i}) \, |v_{i}\rangle\langle v_{i}| & d \text{ even,} \end{cases}\end{split}\]

in the block where the signal qubit, which is the first qubit of the circuit, and the auxiliary qubits start and end in \(|0\rangle\) [3]. The circuit queries \(U\) and \(U^{\dagger}\) alternately \(d\) times and, unlike QSP, never controls them. Its second argument, nancillas, is the number \(a\) of auxiliary qubits that block_encoding_circuit returns.

qsvt_phases takes the coefficients \(c_{k}\) of the Chebyshev series \(P(x) = \sum_{k} c_{k} T_{k}(x)\), with \(T_{k}\) the Chebyshev polynomial of the first kind and degree \(k\). As an example, we use \(P(x) = (3x - x^{3}) / 2\). Since \(x = T_{1}(x)\) and \(x^{3} = \left( 3 T_{1}(x) + T_{3}(x) \right) / 4\), its coefficients are \((0, 9/8, 0, -1/8)\). The polynomial is odd, satisfies \(|P(x)| \le 1\) and \(P(1) = 1\), and amplifies small singular values by a factor \(3/2\).

We apply it to a complex matrix that is not Hermitian, with both block encodings.

[ ]:

matrix = backend.cast(
    [
        [0.30 + 0.10j, 0.20, -0.10j, 0.05],
        [0.10, -0.25 + 0.20j, 0.15, 0.00],
        [0.00, 0.10j, 0.35, -0.20],
        [0.15 - 0.05j, 0.00, 0.10, 0.20j],
    ],
    dtype="complex128",
)
phases = qsvt_phases([0.0, 9 / 8, 0.0, -1 / 8])

for method in ("dense", "lcu"):
    block_encoding, alpha, nauxiliary = block_encoding_circuit(matrix, method=method)
    qsvt = qsvt_circuit(block_encoding, phases, nauxiliary)

    left, values, right = backend.singular_value_decomposition(matrix / alpha)
    transformed = (3 * values - values**3) / 2
    target = backend.matmul(left * transformed, right)

    block = qsvt.unitary(backend)[:4, :4]
    error = backend.matrix_norm(block - target, order=2)
    print(
        f"{method}: {qsvt.nqubits} qubits, {nauxiliary} auxiliary, alpha = {alpha:.3f}, "
        f"error below 1e-7: {error < 1e-7}"
    )
    print("  singular values:", [round(float(value), 3) for value in values])
    print("  transformed:    ", [round(float(value), 3) for value in transformed])

dense: 4 qubits, 1 auxiliary, alpha = 1.000, error below 1e-7: True
  singular values: [0.54, 0.364, 0.295, 0.155]
  transformed:     [0.731, 0.522, 0.429, 0.231]
lcu: 7 qubits, 4 auxiliary, alpha = 1.151, error below 1e-7: True
  singular values: [0.469, 0.316, 0.256, 0.135]
  transformed:     [0.652, 0.459, 0.376, 0.201]

Both circuits implement the same polynomial of the singular values of \(A / \alpha\). The singular values differ between the methods because the normalizations do: the dense method uses \(\alpha = 1\), since the matrix has a norm below one, and the LCU method uses the sum of the magnitudes of the Pauli coefficients.

4. QSP and QSVT side by side#

	QSP (`qsp_circuit`)	QSVT (`qsvt_circuit`)
Input	walk \(W\), built from a Hermitian block encoding	any block encoding \(U\)
Acts on	eigenphases \(\theta\) of \(W\), with \(\sin\theta = \lambda / \alpha\)	singular values \(\varsigma\) of \(A / \alpha\)
Target	complex trigonometric series \(A(\theta) + i C(\theta)\)	real polynomial with definite parity, \(\lvert P \rvert \le 1\)
Queries	\(N\) walk queries, all controlled	\(d\) queries of \(U\) and \(U^{\dagger}\), not controlled
Extra qubit	the QSP qubit, first	the signal qubit, first
Typical use	time evolution \(e^{-i t H}\)	amplification, inversion, filtering of singular values

Practical notes#

Normalization. The circuits process \(A / \alpha\), not \(A\). For the time evolution use \(\tau = \alpha t\), and for QSVT choose the polynomial in terms of \(\varsigma / \alpha\).
Qubit order. Auxiliary qubits come first. The QSP and QSVT circuits add one more qubit before them, so the block of interest is the one where all of those qubits are in \(|0\rangle\).
Exact gates. The walk of QSP is applied under control, so it has to be exact including its global phase. This is why \(-i\) is a gate in qubitized_walk, and why the circuits \(U_{i}\) of lcu_circuit must be exact too.
Precision of the phases. For \(N \gtrsim 30\), qsp_phases loses precision with its default extraction="phase_sums". Use extraction="layer_stripping" for larger degrees.
Cost. block_encoding_circuit needs the full matrix, so it is meant for small matrices and for validation. For larger problems, build the block encoding directly with lcu_circuit.

References#

1. 1. Low and I. L. Chuang, Optimal Hamiltonian Simulation by Quantum Signal Processing, Phys. Rev. Lett. 118, 010501 (2017).
1. 1. Low and I. L. Chuang, Hamiltonian Simulation by Qubitization, Quantum 3, 163 (2019).
1. Gilyén, Y. Su, G. H. Low, and N. Wiebe, Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics, Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing (STOC ’19), 193-204 (2019).
1. 1. Childs and N. Wiebe, Hamiltonian simulation using linear combinations of unitary operations, Quantum Information & Computation 12, 901-924 (2012).
1. 1. Martyn, Z. M. Rossi, A. K. Tan, and I. L. Chuang, Grand Unification of Quantum Algorithms, PRX Quantum 2, 040203 (2021).