Function Differentiation#

The approximation of differential operators is key to efficiently solve partial differential equations. We provide three mechanisms to do it, which are described below.

finite_differences_mpo(order, interval[, ...])

Finite differences operator with noise resilience.

fourier_derivative_mpo(order, interval[, ...])

Nth-order derivative operator via Quantum Fourier Transform.

hdaf_derivative_mpo(order, interval[, M, ...])

Constructs a Matrix Product Operator (MPO) of Hermite Distributed Approximating Functionals (HDAFs).

Finite Differences#

The second-order finite difference method approximates the first and second derivatives of \(f(x)\) as

\[\frac{\partial f(x)}{\partial x} = \frac{f(x + \Delta x) - f(x - \Delta x)}{2 \Delta x} + O(\Delta x^2),\]
\[\frac{\partial^2 f(x)}{\partial x^2} = \frac{f(x + \Delta x) - 2 f(x) + f(x - \Delta x)}{ \Delta x^2} + O(\Delta x^2).\]

MPO representation#

This is translated into the quantum register representation using the displacement operator

\[\begin{split}\hat{\Sigma}^+|s\rangle = \left\{ \begin{array}{ll} |s+1\rangle, & s < 2^{n} \\ 0 & \mbox{else.} \end{array} \right. \quad \hat{\Sigma}^- = \left(\hat{\Sigma}^+\right)^\dagger,\end{split}\]

leading to

\[\ket{\partial_{x}f^{(n)}} \simeq \frac{1}{2\Delta{x}}\left(\hat{\Sigma}^+-\hat{\Sigma}^-\right)\ket{f^{(n)}},\]
\[\ket{\partial^2_{x}f^{(n)}} \simeq \frac{1}{\Delta{x}^2}\left(\hat{\Sigma}^+-2\mathbb{I}+\hat{\Sigma}^-\right)\ket{f^{(n)}}.\]

These MPOs have a bond dimension equal to the number of displacements (\(\chi=3\) for the second derivative), and an error that decays algebraically with the grid step \(O((\Delta x)^2)\) and exponentially in the number of qubits \(O(2^{-2n})\).

Higher-order stencils#

This approximation is improved with smoother formulas to avoid noise artifacts, following Holoborodko. Higher-order derivatives and higher-order stencils can be constructed using the same formalism, with the order parameter controlling the accuracy of the approximation.

An example on how to use these functions is shown in Differentiation.ipynb.

Hermite Distributed Approximating Functionals (HDAFs)#

HDAFs are well-tempered approximations of the Dirac-delta distribution constructed with Hermite polynomials [HNSK91]. They allow approximating operators that are functions of derivatives, making them particularly useful for computing differential operators and time evolution propagators.

Mathematical formulation#

Given a function \(f(x)\), the HDAF approximation reconstructs it through convolution with a smoothed delta function:

\[f(x) \simeq \int \delta_M(x - x';\sigma)f(x')\mathrm{d}x'\]

where the HDAF kernel is:

\[\delta_M(x;\sigma) = \frac{\exp\left(\frac{-x^2}{2\sigma^2}\right)}{\sqrt{2\pi}\sigma} \sum_{m=0}^{M/2}\left(-\frac{1}{4}\right)^m\frac{1}{m!} H_{2m}\left(\frac{x}{\sqrt{2}\sigma}\right)\]

and \(H_n(x)\) is the \(n\)-th Hermite polynomial.

Derivative approximation#

For a derivative of order \(k\), the HDAF approximation becomes:

\[|\partial_x^k f\rangle \approx \left[\Delta x \delta_{M}^{(k)}(0; \sigma) \mathbb{I} + \sum_{i=1}^{W} \Delta x \delta_{M}^{(k)}(i\Delta x; \sigma) \left(\hat{\Sigma}^{+i} + (-1)^{k}\hat{\Sigma}^{-i}\right)\right] |f\rangle,\]

where \(W\) is a suitable cutoff index where the summand vanishes, and

\[\delta_{M}^{(k)}(x; \sigma) = \left(\frac{-1}{\sqrt{2}\sigma}\right)^{k} \frac{\exp\left(\frac{-x^{2}}{2\sigma^{2}}\right)}{\sqrt{2\pi}\sigma} \sum_{m=0}^{M/2}\left(-\frac{1}{4}\right)^{m} \frac{H_{2m + k}\left(\frac{x}{\sqrt{2}\sigma}\right)}{m!}\]

The vectors of HDAF coefficients decay rapidly, so this expansion can be truncated to finite order, yielding an efficient linear combination of displacement operators.

Kinetic propagator#

Beyond differentiation, HDAFs can also approximate the kinetic propagator \(\exp(-it\partial_x^2/2)\) directly in the coordinate basis. This is particularly useful for split-step time evolution methods (see split_step()), where the kinetic propagator is represented as a banded MPO constructed from discrete shift operators, with controllable accuracy determined by the truncation order \(M\).

hdaf_mpo(num_qubits, dx[, M, s0, time, ...])

Constructs a Matrix Product Operator (MPO) of Hermite Distributed Approximating Functionals (HDAFs).

Fourier approximation#

Spectral differentiation based on Fourier techniques achieves significantly higher accuracy than finite differences and allows computing any analytical function \(G(\partial_x)\) of the differential operator.

Mathematical formulation#

Ref. García-Ripoll [GR21] shows that the quantum Fourier transform \(\mathcal{F}\) can be used to construct a differential operator as:

\[G(\partial_x) \to \mathcal{F}^{-1} G(i\mathbf{p}) \mathcal{F}\]

where \(\mathbf{p}\) labels the MPS encodings of the vector of frequencies in Fourier space.

Convergence properties#

For analytic and periodic functions, spectral differentiation achieves exponential convergence \(O(e^{-rN})\) with the number of grid points \(N=2^n\). This implies a super-exponential convergence with the number of qubits: \(O(e^{-r2^n})\).

However, when the function is not periodic, or the interval size does not match the period, the Fourier expansion may lead to localized Gibbs oscillations that decay only algebraically with the number of qubits: \(O(N^{-\alpha}) \sim O(2^{-n\alpha})\).

Implementation notes#

While the momentum operator \(\hat{p}\) admits an MPO representation with bond dimension \(\chi=2\), a naive MPO construction of the QFT leads to a bond dimension that grows linearly with the number of sites. Efficient implementations exploiting bit-reversal symmetry significantly reduce this cost [CSW23].

This can be obtained combining the p_to_n_mpo() function with SeeMPS’s QFT (see Quantum Fourier Transform). The resulting MPOList object may need to be simplified using simplify_mpo() for efficient application.

See also#