Direct interpolative construction of the discrete Fourier transform as a matrix product operator

The quantum Fourier transform (QFT) [2] is one of the most important algorithmic primitives in quantum computing. It is often viewed as providing exponential speedup compared to the fast Fourier transform (FFT), due to its gate complexity of $O(n^{2})$ acting on an $n$-qubit quantum state. Such a state can be viewed as an $N$-dimensional vector with $N=2^{n}$, and under this identification the QFT corresponds to a discrete Fourier transform (DFT), which costs $O(N\log N)=O(2^{n}n)$ operations via the FFT.

However, if one makes use of low-rank structure of the input state, the QFT can often be simulated classically using efficient operations with matrix product states (MPS) and matrix product operators (MPO) [3, 4, 5, 6, 7]. Equivalently, the DFT of functions compressible as quantized tensor trains (QTT) [8, 9] can be computed with exponential speedup. For example, one can remove the long-range quantum gates in the QFT circuit that are close to identity, resulting in the approximate QFT (AQFT) [10]. The AQFT automatically admits an MPO representation due to the bounded interaction length of the circuit [11, 12, 13], and therefore it can be simulated classically on any input states represented as an MPS. The full QFT circuit can also be viewed as the composition of $n$ MPOs of rank 2, which can be applied sequentially to MPS inputs with intermediate compressions [14, 15]. This approach has been called the “superfast Fourier transform” or QTT-FFT [15]. Recently, the entire QFT operator was proven to be compressible as an MPO [1]. In practice, the MPO representation of the QFT, referred to as the QFT MPO, can be obtained by contracting and compressing the entire QFT circuit. Interestingly, for a given error tolerance, the QFT MPO was observed to have smaller ranks compared to the AQFT MPO [1, 16], implying that the long-range gates actually play a role in reducing the entanglement.

Classical simulation of the QFT can be used as a classical algorithm that outperforms the FFT in restricted settings. Concretely, one converts an $N$-dimensional vector into an $n$-site MPS and contracts it with the QFT circuit or the QFT MPO. In this approach, even if the state admits compression as an MPS, the conversion dominates the time complexity if the state is formed as a dense vector and compressed via the TT singular value decomposition (TT-SVD) [8]. A randomized variant of TT-SVD yields an $O(\log N)$ speed-up [1], and other heuristic approaches such as TT-cross [17, 18] can yield $O(N)$ speedup, enabling the QTT-FFT of [15]. If the input vector is the quantized representation of a function that is well approximated in a multiresolution polynomial basis, a direct analytical construction based on the interpolative decomposition can be applied with quantitative control over the error [19]. The interpolating cores introduced in [19] are an essential ingredient in this work. In some situations, the input is automatically available in low-rank tensor format, such as the evolving state of a differential equation solved within the QTT format [20, 21, 22, 23] or a quantum simulation solved using an MPS representation of the state [24, 25, 26].

In these applications, an efficient and accurate construction of the QFT MPO is crucial. While it can be constructed efficiently by contracting the QFT circuit, this approach has a few drawbacks. First, to avoid exponential scaling in $n$, one must perform intermediate compressions before the final MPO is formed, and no rigorous a priori estimate of the accumulated error is available, even though the approach is highly accurate in practice and heuristic arguments support its accuracy. Second, the circuit contraction takes $O(r^{3}n^{2})$ time, where $r$ is the bond dimension, while the MPO ultimately contains at most $O(r^{2}n)$ information, so the computational complexity does not appear to be optimal, as we confirm in this work. Third, the parameters in the circuit decay like $O(1/N)$, so for large $N=2^{n}$ we unavoidably lose accuracy in the long-range gates due to the limitations of numerical precision. This loss of accuracy resembles the approximation imposed by the AQFT circuit, which increases the bond-dimension by the aforementioned observation in [1, 16].

In this paper, we present a simple closed-form construction of the QFT MPO, avoiding circuit contraction and eliminating all three drawbacks mentioned above. The construction involves an interpolative decomposition based on polynomial interpolation on the Chebyshev-Lobatto grid, applied recursively on each qubit, and yields a nearly optimal MPO compression for a given bond dimension. The bond dimension of the construction corresponds to the number of interpolation points. For a fixed number of qubits, the error decreases super-exponentially with respect to the bond dimension, matching the bounds for the Schmidt coefficients derived in [1]. Furthermore, the tensor cores are all identical (excepting the first and last). Finally, we show that the AQFT circuit can also be interpreted as an MPO constructed in the same fashion but with a different interpolation scheme, i.e., piecewise constant interpolation. This interpolation scheme converges more slowly with respect to the number of interpolation points, explaining the observation that the AQFT has a higher MPO rank than the QFT MPO, for a given error tolerance.

We outline the paper as follows. In Section 2 we present some preliminary definitions and notation. In Section 3 we demonstrate how the interpolative decomposition yields a rank bound for the unfolding matrices of the QFT, similar to that of [1]. (In fact, our rank bound is not implied by [1] since it is defined in terms of the infinity norm error, rather than the Frobenius norm error.) We also describe the connection of the rank bound for the unfolding matrices with the complementary low-rank structure [27] that has been observed for the DFT and other operators. This section can also be viewed as a warm-up for Section 4, in which we present our construction of the QFT as an MPO using the interpolative decomposition. In Section 5, we bound the error of this construction. In Section 6, we illustrate the connection to the AQFT.

Consider discrete indices $s,t\in\{0,1,\ldots,N-1\}$, where $N=2^{n}$, and define the matrix of the discrete Fourier transform (DFT)

We comment that in the quantum computing literature, the convention for the QFT is $\exp(2\pi ist/N)$ rather than $\exp(-2\pi ist/N)$, but we use the latter to match the standard definition of the DFT.

We can place $s$ and $t$, respectively, in bijection with bit strings $\sigma_{1:n}=(\sigma_{1},\ldots,\sigma_{n})$ and $\tau_{1:n}=(\tau_{1},\ldots,\tau_{n})$ in $\{0,1\}^{n}$ via

These bits can be viewed as binary digits for $s$ and $t$ with opposite orderings.

Then we can view $\mathbf{F}$ as a tensor $F=F(\sigma_{1:n},\tau_{1:n})$ of size $2^{2n}$ via this identification:

In the QFT, $\tau_{k}$ and $\sigma_{k}$ correspond to the input and output indices of the $k$-th qubit respectively. Notice that $\tau_{1}$ is the least significant bit while $\sigma_{1}$ is the most significant bit, and it is a key point in [1] that only this ordering induces low-rank structure of the QFT as a matrix product operator or quantized tensor train operator, which is defined as the approximate factorization

where we used the notation $[r_{k}]=\{0,\dots,r_{k}-1\}$. The integers $r_{k}$ are often called the bond dimensions or TT ranks, and for simplicity we assume that $r_{k}=r$ for some fixed $r$. The tensors $A_{k}^{\alpha_{k-1},\alpha_{k}}(\sigma_{k},\tau_{k})$ are referred to as tensor cores.

We also define the $m$-th unfolding matrix of $F$ as the $2^{2m}\times 2^{2(n-m)}$ matrix

where we view the multi-indices before and after the semicolon, respectively, as row and column indices of a matrix. It is shown in [8] that a QTT with good approximation exists if the $m$-th unfolding matrix is numerically low-rank for all $m$.

To measure the error of our decompositions it is useful to define a tensor infinity norm

as the largest magnitude of any element of the tensor. We will use this notation for tensors of different shapes, but ultimately the accuracy of our MPO approximation itself is most naturally expressed in this norm.

We also define notations for our interpolation scheme. We will denote the Chebyshev-Lobatto grid points (shifted and scaled to the interval $[0,1]$) by

for $\alpha=0,...,K$. Throughout, we use $K$ to denote grid size (or more properly, the grid size minus one). The corresponding Chebyshev-Lobatto interpolation of a function $f(x)$, i.e., the Lagrange polynomial interpolation on the Chebyshev-Lobatto grid, is

where the $P^{\alpha}$ are degree-$K$ polynomials such that $P^{\alpha}(c^{\beta})=\delta_{\alpha,\beta}$, i.e., the Lagrange interpolating polynomials for this grid. In the case of Chebyshev-Lobatto interpolation, the $P^{\alpha}$ are sometimes called the Chebyshev cardinal functions, which can be evaluated directly with simple trigonometric formulas [28], bypassing the application of the more general formula for Lagrange basis functions.

Although the rank of $F(\sigma_{1:n},\tau_{1:n})$ as an MPO was already bounded in [1], we provide an alternative rank bound from an interpolative perspective. In fact, this point of view yields entrywise control over the error of $F$, which is not captured by the bounds on the Schmidt spectrum in [1].

Our motivation arises from computing

Now if $k\leq m$ and $l>m$, then $2^{-k}2^{l}\sigma_{k}\tau_{l}\in 2\mathbb{Z}$, and therefore

It follows that we can write

where we have overloaded the notation for $F$ slightly to denote the quantized tensor representations of DFT matrices of arbitrary sizes.

It is then useful to define

so

Then we are motivated to consider, for arbitrary $y\in[0,1]$, the function

We can approximate $f_{x}$ with Chebyshev-Lobatto interpolation to obtain

If this equality holds with pointwise error bounded by $\varepsilon$, then the rank-$(K+1)$ approximation of the $m$-th unfolding matrix of $F$ has entrywise error bounded by $\varepsilon$.

The above arguments are reviewed diagrammatically in Fig. 3.1. Next we state and prove an error bound.

The preceding argument does not directly furnish an approximate presentation of $F$ as an MPO. Now we turn to the task of such a construction, which is illustrated diagrammatically in Fig. 4.1.

Recall our definition (3.1) for $y_{\leq m}$ in terms of the bit substring $\tau_{1:m}$, which we reproduce here as

for every $m=1,\ldots,n$. Note that $y_{\leq k}\in[0,1]$ for all $k$, and moreover we can construct the sequence $y_{\leq m}$ recursively via the relation

which holds for $m=1,\ldots,n-1$. Note that $y_{\leq 1}=\frac{1}{2}\tau_{1}$.

We will construct inductively a sequence of tensors $F_{m}\in\mathbb{C}^{2^{m}\times 2^{m}\times(K+1)}$, $m=1,\ldots,n$, with entries specified by

where we point out that $y_{\leq m}$ is always understood to be implicitly defined in terms of $\tau_{1:m}$.

For the base case, the tensor $F_{1}\in\mathbb{C}^{2\times 2\times(K+1)}$ can be written explicitly

Alternatively, we may view $F_{1}^{\beta}$ for each $\beta$ as the $2\times 2$ matrix

Moreover, we can obtain our target $F(\sigma_{1:n},\tau_{1:n})$ as

under the convention that $\alpha=0$ is the index of the leftmost node $c^{\alpha}=0$.

Then we claim that $F_{m+1}$ can be obtained approximately from $F_{m}$ by attaching a single tensor core on the right. To see this, expand:

Now for any $k=1,\ldots,m$, observe that $2^{-k}2^{m+1}\sigma_{k}\in 2\mathbb{Z}$, so the second factor in the last expression is just $1$. Moreover, we can substitute $\sum_{l=1}^{m}2^{l-m-1}\tau_{l}=y_{\leq m}$ in the third factor and $y_{\leq m+1}=\frac{1}{2}y_{\leq m}+\frac{1}{2}\tau_{m+1}$ in the last factor to obtain

where we retain our definition of $f_{y}:[0,1]\rightarrow\mathbb{R}$ from the last section.

Then we can insert an interpolative decomposition

into our expression for $F_{m+1}^{\beta}(\sigma_{1:m+1},\tau_{1:m+1})$, yielding

where we have defined the tensor core $A\in\mathbb{C}^{2\times 2\times(K+1)\times(K+1)}$ via

Alternatively we can write

Then we have argued that we can approximate $F$ as an MPO with leftmost core $A_{\mathrm{L}}=F_{1}$, internal cores $A$, and rightmost core $A_{\mathrm{R}}\in\mathbb{R}^{2\times 2\times(K+1)}$ given by

Since $\sum_{\alpha=0}^{K}P^{\alpha}(x)=1$ for all $x\in[0,1]$, we can neatly express the leftmost core in terms of $A$ as

The following result bounds the entrywise error of the MPO constructed as above.

The AQFT [10] on $n$ qubits with approximation level $b\in\{0,\ldots,n-1\}$ is defined by

up to the choice of convention for normalization and complex conjugation. Relative to the formula (2.2) for the QFT/DFT, in the AQFT we discard all terms in the exponent with $k-l>b$. If $b=n-1$ we recover the exact QFT, and if $b=0$ we get the Hadamard transform. The quantum circuits for the QFT and AQFT are shown in Fig. 6.2. We prove the following theorem regarding the connection to interpolative decomposition: