1 Introduction↩︎

Quantum computers are expected to solve certain computational problems more efficiently than classical computers by exploiting high-dimensional quantum state spaces and coherent quantum dynamics. Potential applications include scientific computing, quantum simulation, optimization, finance, and machine learning; for instance, see the surveys [1], [2].

Many quantum algorithms for linear algebra and differential equations can be understood as eigenvalue transformations. If a matrix \(A\in\mathbb{C}^{N\times N}\) is diagonalizable as \(A=VDV^{-1}\) and \(f\) is defined on the spectrum of \(A\), the eigenvalue transformation associated with \(f\) is \[f(A):=Vf(D)V^{-1}.\] For analytic \(f\), this definition agrees with the usual matrix functional calculus and extends to non-diagonalizable matrices by continuity. The corresponding quantum task is to block-encode \(f(A)\) in a quantum circuit, or to prepare a normalized state proportional to \(f(A)\ket{\psi}\) for an input state \(\ket{\psi}\). The well-known QSVT framework [3] solves this problem for Hermitian \(A\), but it does not apply to general matrices.

For the important special case \(f(z)=e^{-tz}\) on the right half-plane \(\{\mathrm{Re}(z)\ge 0\}\), the linear combination of Hamiltonian simulation (LCHS) approach is particularly effective. For a matrix \(A\), write its Cartesian decomposition as \(A=X+\mathrm{i}Y\), where \(X\) and \(Y\) are Hermitian. Under the dissipativity assumption \(X = \frac{A+A^\dagger}{2}\succeq 0\), prior works [4]–[6] prove formulas of the form \[\label{eq:32first32lchs32formula} e^{-tA} = (2\pi)^{-1/2}\int_{-\infty}^\infty \hat{h}(u) e^{-\mathrm{i}t(uX+Y)} \,\mathrm{d}u\tag{1}\] for suitable kernels \(\hat{h}\). Explicit examples include \[\label{eq:32LCHS32exp40-tA41} \hat{h}(u) = \sqrt{\frac{2}{\pi}}\frac{1}{1+u^2}, \quad \text{or}\quad \hat{h}(u) = \frac{e^{2^\beta-(1+\mathrm{i}u)^\beta}}{\sqrt{2\pi}(1-\mathrm{i}u)},\quad 0<\beta<1.\tag{2}\] These identities are nontrivial generalizations of the corresponding scalar identities, because \(X\) and \(Y\) generally do not commute. On a quantum computer, the unitary factors \(e^{-\mathrm{i}t(uX+Y)}\) can be implemented by well-studied Hamiltonian simulation algorithms [3], [7], [8]. After truncation and discretization, 1 yields quantum circuits that encode \(e^{-tA}\) to additive precision \(\epsilon\) using \(\mathcal{O}\left(t\log^{1+\mathcal{O}\left(1\right)}\frac{1}{\epsilon}\right)\) queries to block-encodings of \(A\) and \(A^\dagger\).

Ref. [9] observed that certain truncated kernels \(\hat{h}\) need not satisfy the exact identity 1 , but can still satisfy an approximation of the form \[\label{eq:32second32lchs32formula} e^{-tA} \approx (2\pi)^{-1/2}\int_{-R}^R \hat{h}(u) e^{-\mathrm{i}t(uX+Y)} \,\mathrm{d}u\tag{3}\] for some finite \(R\). This relaxation leads to quantum algorithms with improved query complexity \(\mathcal{O}\left(t\log\frac{1}{\epsilon}\right)\).

From the viewpoint of [6], [9], the above LCHS formulas can be interpreted as follows: one first finds a function \(h(x)\) that approximates \(e^{-x}\) on \([0,\infty)\), and then uses the Fourier transform \(\hat{h}\) as the LCHS kernel. A natural question is whether this perspective can be generalized to other functions \(f(z)\).

For \(e^{-z}\), the scalar identity \(e^{-(x+\mathrm{i}y)}=e^{-x}e^{-\mathrm{i}y}\) separates the decay in \(x\) from the phase in \(y\), and the phase factor is incorporated into the Hamiltonian simulation term. For a general analytic function, such a separation is unavailable. It is therefore natural to Fourier transform in both the real and imaginary directions, leading to the following question:

If \(h(x,y)\) approximates \(f(x+\mathrm{i}y)\) in some region, will the formula \[\label{eq:32lchs32formula} f(A) \approx (2\pi)^{-1}\int_{u^2+v^2\le R^2} \hat{h}(u,v) e^{-\mathrm{i}(uX+vY)} \,\mathrm{d}u \,\mathrm{d}v\qquad{(1)}\] hold in some sense?

In this paper, we give a positive answer to this question using Weyl calculus.

Even with this representation, the continuous LCHS formula ?? is not yet directly suitable for coherent quantum circuit implementation. In practice, one usually needs a discrete LCHS formula \[\label{eq:32disc32lchs32formula} f(A) \approx \sum_{(u_j,v_j)\in \mathcal{G}\subset D_R} w_{j}\, e^{-\mathrm{i}(u_jX+v_jY)}\tag{4}\] where \(\mathcal{G}\) is a grid contained in the disc \(D_R=\{(u,v) \mid u^2+v^2 \le R^2\}\). In previous LCHS works, such discrete formulas are obtained by truncating the continuous formula and then applying quadrature rules. This approach has two limitations. First, explicitly constructed kernels \(\hat{h}\) may have suboptimal constants or scalings, and this inefficiency is inherited by the discretized formula. Second, quadrature rules are most effective for smooth or analytic kernels. If \(\hat{h}\) is not sufficiently smooth, even when \(\hat{h}\) decays rapidly, the number of discretized terms can become large. Based on our Weyl calculus framework, we develop a convex optimization method that directly finds good discrete formulas of the form 4 under a prescribed resource budget. Interestingly, applying this optimizer to \(f(z)=e^{-z}\) suggests a closed-form kernel for 3 that improves the cost relative to prior art [9].

For the quantum implementation of 4 , the cost is largely determined by two parameters. The first is the truncation radius \(R\), which controls the maximum simulation time in the Hamiltonian simulation circuits. The second is the normalization factor \(\lambda = \sum_j |w_j|\) (or \(\lambda = (2\pi)^{-1} \mathchoice {\left\lVert\hat{h}\right\rVert} {\lVert\hat{h}\rVert} {\lVert\hat{h}\rVert} {\lVert\hat{h}\rVert} _{L^1(\mathbb{R}^2)}\) in the continuous case ?? ), which controls the postselection probability in LCU-based implementations. Amplitude amplification [10] can compensate for this normalization but at the cost of proportionally increasing the simulation time. Therefore, this paper uses \[\Lambda = \lambda R\] as the cost function of an LCHS formula.

1.1 Main results↩︎

1.1.1 Weyl-calculus framework for general LCHS↩︎

We develop a general LCHS theory based on Weyl calculus, providing a systematic way to construct LCHS formulas for analytic matrix functions.

Moreover, the framework allows \(\hat{h}\) to be a finite measure, such as a Dirac comb, so it directly covers discrete weight distributions.

1.1.2 Closed-form LCHS constructions for analytic functions↩︎

Using the Weyl calculus framework, we construct LCHS formulas for analytic functions and derive explicit bounds on the cost function \(\Lambda\).

We also provide a sharper upper bound in 7 by replacing \(S\) with the local Wiener norm \(\mathchoice {\left\lVert f\right\rVert} {\lVert f\rVert} {\lVert f\rVert} {\lVert f\rVert} _{W(\mathbb{D}_{1+\rho})}\), defined in 7 . However, this sharper bound is less explicit because the local Wiener norm is usually difficult to compute exactly. Additionally, if the assumption \(\mathchoice {\left\lVert A\right\rVert} {\lVert A\rVert} {\lVert A\rVert} {\lVert A\rVert} \le 1\) is replaced by a more general numerical range assumption, we give a similar construction in 9.

1.1.3 Improved prefactor for dissipative dynamics↩︎

For the specific function \(f(z)=e^{-z}\) (\(\mathrm{Re}(z)\ge 0\)), we find a kernel function that improves the cost by a factor of \(2.1\times\) over prior art [9] and is within a factor of \(1.5\times\) of the theoretical lower bound. As in previous LCHS works, our construction also generalizes to the time-dependent ODE \[\frac{\,\mathrm{d}}{\,\mathrm{d}t} \vec{\psi}(t) = -A(t)\vec{\psi}(t),\] which is solved by the time-ordered operator \[\label{eq:32time32ordered32operator32def} \mathcal{T}e^{-\int_0^t A(s) \,\mathrm{d}s}: \vec{\psi}(0) \mapsto \vec{\psi}(t).\tag{5}\]

1.1.4 Direct optimization of discrete LCHS weights↩︎

We develop a convex optimization framework for finding good discrete weights \(\{w_j\}\) in 4 under a given budget \((R,\lambda)\). The theoretical guarantee is the discrete analogue of 1, corresponding to the case where \(h\) is a finite sum of Fourier modes. The resulting programs are relatively tight because they optimize directly over the weights \(\{w_j\}\), rather than over a parametrized ansatz for the weight coefficients. The costs of the optimized LCHS formulas are typically much smaller than those of the closed-form construction in 2. For a fixed approximation region \(\Omega\), we use a fixed-stepsize grid \(\mathcal{G}\); therefore, the number of terms in 4 only grows like \(\mathcal{O}\left(R^2\right) = \mathcal{O}\left(\log^2\frac{1}{\epsilon}\right)\).

For the one-dimensional interval \([0,T]\) and the two-dimensional regions \([-1,1]^2\) and \(D_1 = \{(x,y):x^2+y^2\le 1\}\), we derive concrete local Wiener norm upper bounds and present the corresponding optimization problems. We test the optimization method on several examples.

1.2 Related work↩︎

The closest line of work is the linear combination of Hamiltonian simulation (LCHS) approach to non-unitary dynamics. The original LCHS method [4] represents \(e^{-tA}\), or the more general time-dependent version 5 , for \(A=L+\mathrm{i}H\) with \(L\succeq 0\), as a linear combination of Hamiltonian evolutions using the formula 2 , and achieves optimal dependence on the input state preparation cost. The improved LCHS constructions [5], [6], [9], [11] use different kernel functions to obtain better query complexity dependence on the target precision \(\epsilon\) for linear non-unitary dynamics. Subsequent works have refined this direction by extending the operator setting [12] and developing problem-specific or implementation-oriented variants [13]–[15].

The Lap-LCHS framework [16] extends LCHS from matrix exponentials to functions that admit a suitable Laplace representation, thereby covering several eigenvalue transformations such as powers of shifted inverses and mass-matrix differential equations. The core idea is that if a function \(f(z)\) has a convergent Laplace transform representation \[f(z) = \int_0^\infty g(t) e^{-tz} \,\mathrm{d}t,\quad \mathrm{Re}(z)\ge 0,\] then one may implement \(f(A)\) as a linear combination of terms \(e^{-tA}\), which can be handled by the original LCHS framework. However, only certain analytic functions that decay on the right half-plane admit convergent Laplace transforms. This excludes many common functions, including polynomials \(p(A)\), logarithms \(\log(A)\), and matrix powers \(A^p\;(p>0)\). In contrast, our work treats general analytic functions through a two-dimensional Fourier/Weyl-calculus representation and also allows the weights to be discrete from the outset.

Another general approach to eigenvalue transformation is based on contour integrals. For analytic \(f\), Cauchy’s formula expresses \(f(A)\) as an integral of resolvents \((zI-A)^{-1}\), and quadrature reduces the task to a linear combination of shifted linear-system solves [17]–[19]. Recent work gives a more systematic contour-integral-based quantum eigenvalue transformation framework and a detailed complexity analysis [20]. This route is quite general for non-Hermitian and non-normal matrices, but its cost depends on the contour geometry and the implementation of shifted inverses. The query complexity typically deteriorates when \(A\) is non-diagonalizable or has a poorly conditioned eigenvector basis [20], an issue that does not arise in LCHS-based methods.

Polynomial and signal-processing approaches form a complementary family. Quantum singular value transformation [3] and related works [21]–[23] provide a powerful framework for applying polynomial transformations to singular values, but singular value transformations do not directly give eigenvalue transformations for general non-normal matrices. Quantum eigenvalue processing [24] overcomes this obstruction for polynomial eigenvalue transformations using Chebyshev or Faber history states. Its complexity also depends on the condition number of the Jordan basis of \(A\), similar to the contour-integral approach. Recent work [25] presents a conceptually simpler eigenvalue-transformation framework based on GQSP [22] and the “regular block encoding” of the matrix \(A\).

1.3 Organization of the paper↩︎

The rest of the paper is organized as follows. 2 fixes the Fourier-transform convention and introduces the Wiener norm used to measure LCHS approximation errors. In 3, we develop the Weyl-calculus framework for two-variable Fourier representations and prove the main LCHS error bound (1). In 4, we use this bound to construct explicit LCHS formulas for general analytic functions (2). In 5, we adapt the Weyl-calculus argument to time-ordered non-unitary dynamics and then give the improved kernel function for dissipative dynamics in 5.1 (3). In 6, we present the discrete-weight optimization framework and describe the corresponding convex programs for a one-dimensional interval, a two-dimensional square, and a two-dimensional disc. 7 summarizes coherent and hybrid quantum implementations and explains how the LCHS parameters \((R,\lambda)\) enter the quantum resource estimates. Finally, 8 discusses potential future directions.

2 Preliminaries: Fourier transforms and Wiener norms↩︎

We use the following convention for the Fourier transform and its inverse: \[g(\mathbf{x})=(2 \pi)^{-n / 2} \int_{\mathbb{R}^n} \hat{g}(\boldsymbol{\xi}) e^{-\mathrm{i}\boldsymbol{\xi}\cdot \mathbf{x}} \,\mathrm{d}\boldsymbol{\xi}\] \[\hat{g}(\boldsymbol{\xi})=(2 \pi)^{-n / 2} \int_{\mathbb{R}^n} g(\mathbf{x}) e^{\mathrm{i}\boldsymbol{\xi}\cdot \mathbf{x}} \,\mathrm{d}\mathbf{x}\] To include discrete LCHS formulas in the same framework, we also use the Fourier–Stieltjes convention. Namely, for a finite complex Borel measure \(\mu\), we write \(\hat{g} = \mu\) if \[g(\mathbf{x})=(2 \pi)^{-n / 2} \int_{\mathbb{R}^n} e^{-\mathrm{i}\boldsymbol{\xi}\cdot \mathbf{x}} \, \mathrm{d} \mu(\boldsymbol{\xi})\] holds. In particular, \(\widehat{e^{-\mathrm{i}\mathbf{a}\cdot \mathbf{x}}} = (2 \pi)^{n / 2} \delta_{\mathbf{a}}\), where \(\delta_{\mathbf{a}}\) is the Dirac measure at \(\mathbf{a}\). We refer the reader to 9.1 for more details on the Fourier–Stieltjes algebra. We use the standard notation \[A(\mathbb{R}^n) := \{ g: \hat{g}\in L^1(\mathbb{R}^n)\}, \quad B(\mathbb{R}^n) := \{g: \hat{g} = \mu \text{ such that } \int_{\mathbb{R}^n} \,\mathrm{d} \mathchoice {\left|\mu\right|} {|\mu|} {|\mu|} {|\mu|} \le\infty\}.\] For \(g\in B(\mathbb{R}^n)\) and \(\mu = \hat{g}\), we define the Wiener norm \[\label{eq:32wiener32norm32def} \mathchoice {\left\lVert g\right\rVert} {\lVert g\rVert} {\lVert g\rVert} {\lVert g\rVert} _W := (2 \pi)^{-n / 2} \int_{\mathbb{R}^n} \,\mathrm{d} \mathchoice {\left|\mu\right|} {|\mu|} {|\mu|} {|\mu|} \;\xlongequal{g\in A(\mathbb{R}^n)} (2 \pi)^{-n / 2} \int_{\mathbb{R}^n} \mathchoice {\left|\hat{g}(\boldsymbol{\xi})\right|} {|\hat{g}(\boldsymbol{\xi})|} {|\hat{g}(\boldsymbol{\xi})|} {|\hat{g}(\boldsymbol{\xi})|} \,\mathrm{d}\boldsymbol{\xi}.\tag{6}\] Under this norm, \(B(\mathbb{R}^n)\) is a Banach space. The space \(A(\mathbb{R}^n)\) is the closed subspace corresponding to absolutely continuous finite measures, since an \(L^1\) Fourier transform \(\hat{g}(\boldsymbol{\xi})\) can be identified with the measure \(\mu = \hat{g}(\boldsymbol{\xi})\,\mathrm{d}\boldsymbol{\xi}\).

As an important example, when \(g\) is a periodic function given by a Fourier series, \(\mathchoice {\left\lVert g\right\rVert} {\lVert g\rVert} {\lVert g\rVert} {\lVert g\rVert} _W\) is the \(\ell^1\) norm of its Fourier coefficients.

Suppose that a smooth function \(g\) is specified only on a closed set \(\Omega \subset \mathbb{R}^n\), in the sense that it extends to a \(C^\infty\) function on a neighborhood of \(\Omega\). We define its local Wiener norm by \[\label{eq:32restricted32wiener32norm32def} \mathchoice {\left\lVert g\right\rVert} {\lVert g\rVert} {\lVert g\rVert} {\lVert g\rVert} _{W(\Omega)} :=\;\inf_{\tilde{g}\in B(\mathbb{R}^n),\;\tilde{g}|_\Omega = g}\; \mathchoice {\left\lVert\tilde{g}\right\rVert} {\lVert\tilde{g}\rVert} {\lVert\tilde{g}\rVert} {\lVert\tilde{g}\rVert} _W\tag{7}\] When \(\Omega\) is compact and convex, the infimum can be taken over more regular extensions without changing its value: \[\label{eq:32regular32extension} \mathchoice {\left\lVert g\right\rVert} {\lVert g\rVert} {\lVert g\rVert} {\lVert g\rVert} _{W(\Omega)} =\;\inf_{\tilde{g}\in B(\mathbb{R}^n),\;\tilde{g}|_\Omega = g}\; \mathchoice {\left\lVert\tilde{g}\right\rVert} {\lVert\tilde{g}\rVert} {\lVert\tilde{g}\rVert} {\lVert\tilde{g}\rVert} _W =\;\inf_{\hat{h}\in C^{\infty}(\mathbb{R}^n),\;h\in A(\mathbb{R}^n),\;h|_\Omega = g}\; \mathchoice {\left\lVert h\right\rVert} {\lVert h\rVert} {\lVert h\rVert} {\lVert h\rVert} _W.\tag{8}\] This equality is proved in 9.2.

The Wiener norm can often be controlled by Sobolev norms. Roughly speaking, more than \((n/2)\)-th weak derivatives are enough to guarantee a finite Wiener norm. In two dimensions, Cauchy–Schwarz gives \[\label{eq:32H232estimate} \begin{align} \mathchoice {\left\lVert g\right\rVert} {\lVert g\rVert} {\lVert g\rVert} {\lVert g\rVert} _{W(\mathbb{R}^2)} &=\frac{1}{2\pi} \iint_{\mathbb{R}^2} \mathchoice {\left|\hat{g}(u,v)\right|} {|\hat{g}(u,v)|} {|\hat{g}(u,v)|} {|\hat{g}(u,v)|} \,\mathrm{d}u \,\mathrm{d}v \\ &\le \frac{1}{2\pi} \left( \iint_{\mathbb{R}^2} \frac{1}{(1+u^2+v^2)^2} \,\mathrm{d}u \,\mathrm{d}v \right)^{1/2} \left( \iint_{\mathbb{R}^2} (1+u^2+v^2)^2 \mathchoice {\left|\hat{g}(u,v)\right|} {|\hat{g}(u,v)|} {|\hat{g}(u,v)|} {|\hat{g}(u,v)|} ^2 \,\mathrm{d}u \,\mathrm{d}v \right)^{1/2}\\ &= \frac{1}{2\sqrt{\pi}} \mathchoice {\left\lVert(1-\Delta)g\right\rVert} {\lVert(1-\Delta)g\rVert} {\lVert(1-\Delta)g\rVert} {\lVert(1-\Delta)g\rVert} _{L^2(\mathbb{R}^2)}. \end{align}\tag{9}\] Thus, every \(g\) in the Sobolev space \(H^2(\mathbb{R}^2)\) has finite Wiener norm.

Exact local Wiener norms are usually difficult to compute, but useful upper bounds can be obtained from extensions. If \(g\) is defined only on a closed set \(\Omega\subset\mathbb{R}^2\), we may extend it to a function in \(H^2(\mathbb{R}^2)\) and then apply 9 . The following standard extension result is sufficient for this purpose.

We will also use a cruder \(W^{2,\infty}\) bound for the disc \(D_a\) of radius \(a\ge 1\). The advantage of this estimate is that it gives explicit constants.

The construction of \(E_ag\) and the proof of the estimate are given in 9.3.

3 Weyl calculus and generalized LCHS↩︎

Weyl calculus provides a systematic way to lift scalar functions of several variables to functions of noncommuting operators. It was originally developed in connection with the quantization of functions of position and momentum operators.

We only need the two-variable setting, which already contains the essential ideas. Given a scalar function \(g(x,y)\) and Hermitian matrices \(X,Y\), we want to define an operator-valued function \(\mathrm{Op}\left(g\right)\left(X,Y\right)\). For the monomial \(p_{m,n}(x,y)=x^m y^n\), a natural definition is \[\label{eq:32weyl32monomial} \mathrm{Op}\left(p_{m,n}\right)\left(X,Y\right) = \frac{1}{\binom{m+n}{m}}\sum_{\sigma} A_{\sigma(1)}A_{\sigma(2)}\cdots A_{\sigma(m+n)},\tag{11}\] where the sum is over all binary strings \(\sigma\) containing \(m\) zeros and \(n\) ones, with \(A_0=X\) and \(A_1=Y\). For example, \(\mathrm{Op}\left(p_{2,1}\right)\left(X,Y\right) = \frac{1}{3}(X^2Y+XYX+YX^2)\). We then define \(\mathrm{Op}\left(\cdot\right)\left(X,Y\right)\) for polynomials and power series by linearity. With this definition, one checks that \[\begin{align} \mathrm{Op}\left((ax+by)^n\right)\left(X,Y\right) & = (aX+bY)^n,\tag{12}\\ \mathrm{Op}\left(e^{ax+by}\right)\left(X,Y\right) & = e^{aX+bY}.\tag{13} \end{align}\] Indeed, 12 follows from the expansion \((ax+by)^n = \sum_{k=0}^n \binom{n}{k} a^kb^{n-k}x^ky^{n-k}\): the operator \(\mathrm{Op}\left(x^ky^{n-k}\right)\left(X,Y\right)\) averages over all products with \(k\) copies of \(X\) and \(n-k\) copies of \(Y\), matching the corresponding terms in \((aX+bY)^n\). The identity 13 then follows by applying 12 term by term to the power series \(e^{ax+by} = \sum_{k=0}^\infty \frac{1}{k!}(ax+by)^k\).

For a function \(g\) with \(\mathchoice {\left\lVert\hat{g}\right\rVert} {\lVert\hat{g}\rVert} {\lVert\hat{g}\rVert} {\lVert\hat{g}\rVert} _{L^1(\mathbb{R}^2)} < \infty\), Weyl [28] suggested the definition \[\label{eq:32weyl32fourier32defi} \mathrm{Op}\left(g\right)\left(X,Y\right) = \frac{1}{2\pi}\iint_{\mathbb{R}^2} \hat{g}(u,v) e^{-\mathrm{i}(uX+vY)} \,\mathrm{d}u \,\mathrm{d}v\tag{14}\] The motivation is that \[\begin{align} \frac{1}{2\pi}\iint_{\mathbb{R}^2} \hat{g}(u,v) e^{-\mathrm{i}(uX+vY)} \,\mathrm{d}u \,\mathrm{d}v &= \frac{1}{2\pi}\iint_{\mathbb{R}^2} \hat{g}(u,v) \mathrm{Op}\left(e^{-\mathrm{i}(ux+vy)}\right)\left(X,Y\right)\,\mathrm{d}u \,\mathrm{d}v\\ &= \mathrm{Op}\left(\frac{1}{2\pi}\iint_{\mathbb{R}^2} \hat{g}(u,v) e^{-\mathrm{i}(ux+vy)} \,\mathrm{d}u \,\mathrm{d}v\right)\left(X,Y\right)\\ &= \mathrm{Op}\left(g\right)\left(X,Y\right),\\ \end{align}\] where the second equality uses the exponential identity 13 .

The relevance of Weyl calculus to LCHS comes from the following observation. Let \(f(z)\) be analytic and set \(g(x,y)=f(x+\mathrm{i}y)\). If \(A=X+\mathrm{i}Y\), then 12 suggests that \(\mathrm{Op}\left(g\right)\left(X,Y\right)=f(X+\mathrm{i}Y)\), since both \(g\) and \(f\) can be expanded in power series. At the same time, 14 represents \(\mathrm{Op}\left(g\right)\left(X,Y\right)\) as a linear combination of unitary matrices \(e^{-\mathrm{i}(uX+vY)}\). This is precisely the LCHS form. The main obstacle is that the Fourier transform of an analytic function is usually not an \(L^1\) function or a finite Borel measure, so 14 cannot be applied directly. The rest of this section develops the framework needed to overcome this obstruction.

We first recall a few functional-analytic notions. The Schwartz space \(\mathcal{S}(\mathbb{R}^2)\) consists of smooth functions whose derivatives of all orders decay faster than any polynomial at infinity. A matrix-valued distribution is a continuous linear functional from the test function space \(C^\infty(\mathbb{R}^2)\) to the matrix space \(\mathrm{Mat}_N(\mathbb{C})\). The support of a matrix-valued distribution \(T\), denoted by \(\mathrm{supp}(T)\), is the complement of the largest open set \(U \subset \mathbb{R}^2\) such that \(T(g) = 0\) for all test functions \(g\) supported in \(U\).

The lemma states that \(\mathrm{Op}\left(\cdot\right)\left(X,Y\right)\) is compactly supported, but it does not identify the support. The next result gives the needed geometric description.

By the Toeplitz–Hausdorff theorem, \(W(A)\) is a convex subset of \(\mathbb{C}\). It contains all eigenvalues of \(A\), but it can be larger than the convex hull of the eigenvalues unless \(A\) is normal. Two simple but important examples are

• \(A+A^\dagger \succeq 0 \Longleftrightarrow W(A)\subset \{z: \mathrm{Re}(z)\ge 0\};\)

• \(\mathchoice {\left\lVert A\right\rVert} {\lVert A\rVert} {\lVert A\rVert} {\lVert A\rVert} \le R \Longrightarrow W(A)\subset \mathbb{D}_R\).

The following proposition gives the connection between eigenvalue transformation and Weyl calculus.

3 shows that \(\mathrm{Op}\left(g\right)\left(X,Y\right)\) is well-defined for \(C^\infty\) functions and depends only on the values of \(g\) on a bounded region. However, explicitly computing Fourier transforms of \(C^\infty\) functions is usually difficult. We therefore relax the smoothness requirement and allow \(g\in B(\mathbb{R}^2)\). Let \[\label{eq:32FS32f} g(x,y) = (2 \pi)^{-1} \int_{\mathbb{R}^2} e^{-\mathrm{i}(ux+vy)} \,\mathrm{d}\mu(u,v)\tag{18}\] for some finite complex Borel measure \(\mu\), and we naturally define \[\label{eq:32def32opg32general} \mathrm{Op}\left(g\right)\left(X,Y\right):=(2 \pi)^{-1} \int_{\mathbb{R}^2} e^{-\mathrm{i}(uX+vY)} \,\mathrm{d}\mu(u,v),\tag{19}\] which is a generalization to 14 .

The preceding lemma shows that, for a function \(g\) defined on \(\Omega\), the operator \(\mathrm{Op}\left(g\right)\left(X,Y\right)\) is well-defined whenever \(g\) admits an extension with finite Wiener norm. We will often approximate \(g\) by another function, and the next estimate controls the resulting operator error.

This theorem is the central estimate in the Weyl-calculus approach to LCHS. It reduces the construction of approximate LCHS formulas for an analytic function \(f\) to a scalar approximation problem: one only needs to approximate the two-variable function \(f(x+\mathrm{i}y)\) on the numerical range of \(A\). Since the Fourier–Stieltjes transform of \(h\) may be any finite Borel measure, the same statement covers both continuous and discrete LCHS formulas.

4 LCHS-based eigenvalue transformation for analytic functions↩︎

In this section, we use the Weyl-calculus error bound to construct explicit LCHS formulas for general analytic functions. The construction depends on a region containing the numerical range of the input matrix \(A\). We first treat the common normalization \(\mathchoice {\left\lVert A\right\rVert} {\lVert A\rVert} {\lVert A\rVert} {\lVert A\rVert} \le 1\), which is especially natural in the block-encoding query model (2). Under this assumption, \(W(A)\subset \mathbb{D}_1\). The more general case \(\mathchoice {\left\lVert A\right\rVert} {\lVert A\rVert} {\lVert A\rVert} {\lVert A\rVert} \le \alpha\) can be reduced by scaling. In addition, we will also state a version for arbitrary bounded convex numerical ranges at the end of the section.

Assume that \(f(z)\) is analytic in a neighborhood of \(\mathbb{D}_1\). Then there exists \(\rho>0\) such that \(f\) is analytic and bounded in the interior of \(\mathbb{D}_a\), where \(a=1+\rho\). Define the two-variable function \(g:\mathbb{R}^2\to \mathbb{C}\) by \(g(x,y)=f(x+\mathrm{i}y)\). Our goal is to construct a smooth approximation to \(g\) using only local data on \(D_a\).

We fix an extension \(g_2\) of \(g|_{D_a}\) whose Fourier transform is an \(L^1\) function: \[\label{eq:32g232L132extension32e95disc} g_2=g\quad \text{on }D_a, \qquad \hat{g}_2\in L^1(\mathbb{R}^2).\tag{21}\] At this stage, the construction works with any such extension \(g_2\). Specific choices of \(g_2\) will be used later in 7 and 8.

Put \[\label{eq:32rho32and32M32e95disc} \rho=a-1,\qquad M= \mathchoice {\left\lVert g_2\right\rVert} {\lVert g_2\rVert} {\lVert g_2\rVert} {\lVert g_2\rVert} _W=(2\pi)^{-1} \mathchoice {\left\lVert\hat{g}_2\right\rVert} {\lVert\hat{g}_2\rVert} {\lVert\hat{g}_2\rVert} {\lVert\hat{g}_2\rVert} _{L^1(\mathbb{R}^2)}.\tag{22}\] Then \[\label{eq:32g232Linfty32M32e95disc} \mathchoice {\left\lVert g_2\right\rVert} {\lVert g_2\rVert} {\lVert g_2\rVert} {\lVert g_2\rVert} _{L^\infty(\mathbb{R}^2)}\le M.\tag{23}\]

We convolve \(g_2\) with the Gaussian kernel \[K_\gamma(x,y)=\frac{\gamma^2}{\pi}e^{-\gamma^2(x^2+y^2)}.\] This gives the analytic function \[\label{eq:32g332e95disc} g_3(x,y) = K_\gamma*g_2(x,y).\tag{24}\] Equivalently, in the Fourier domain, \[\label{eq:32g332Fourier32e95disc} \hat{g}_3(\xi,\eta)=e^{-\frac{\xi^2+\eta^2}{4\gamma^2}}\hat{g}_2(\xi,\eta).\tag{25}\] Given an accuracy target \(\epsilon\), we will choose \(\gamma\) and \(R\) so that \(g_3\) satisfies the following three properties:

• \(\mathchoice {\left\lVert g_3\right\rVert} {\lVert g_3\rVert} {\lVert g_3\rVert} {\lVert g_3\rVert} _{W}\) is bounded uniformly for \(\epsilon\).

• \(\mathchoice {\left\lVert g_3-g\right\rVert} {\lVert g_3-g\rVert} {\lVert g_3-g\rVert} {\lVert g_3-g\rVert} _{W(D_1)} \le \frac{\epsilon}{2}\).

• \((2\pi)^{-1} \mathchoice {\left\lVert\hat{g_3}\right\rVert} {\lVert\hat{g_3}\rVert} {\lVert\hat{g_3}\rVert} {\lVert\hat{g_3}\rVert} _{L^1(\mathbb{R}^2\setminus D_R)}\le\frac{\epsilon}{2}\) and \(R = \mathcal{O}\left(\log\frac{1}{\epsilon}\right)\).

Together, these properties allow us to apply 6 with the truncated kernel \(\hat{h}=\hat{g}_3|_{D_R}\), yielding an LCHS formula with controlled approximation error and cost. The next three lemmas make the choices of \(\gamma\) and \(R\) precise.

The preceding theorem is not fully explicit because the extension \(g_2\) is not explicitly constructed. In principle, one can compute such an extension numerically, since finding a minimum-Wiener-norm extension is a convex optimization problem. An alternative choice is to use the explicit extension \(g_2=E_ag\) from 2. This gives a cost bound in terms of derivative bounds for \(f\) up to second order, which we state in the following theorem.

[thm: construction 1 for disc,thm: construction 2 for disc] were stated for the case where the numerical range of \(A=X+\mathrm{i}Y\) is contained in the unit disc \(D_1\). The same construction extends to other numerical ranges. Let \(\Omega\) be the numerical range of \(A\), and let \(\Omega_\rho\) be the \(\rho\)-neighborhood of \(\Omega\). If \(f\) is analytic in \(\Omega_\rho\), we can again choose an extension \(g_2\) with finite Wiener norm and convolve it with a Gaussian kernel. The proofs of 6 and 8 are unchanged. In 7, the only geometric property of \(D_1\) used is that every ball of radius \(\rho\) centered in \(D_1\) is contained in \(D_{1+\rho}\). The corresponding property holds for any compact set \(\Omega\) and its \(\rho\)-neighborhood \(\Omega_\rho\). Finally, the use of 2 can be replaced by the extension result 1, which only changes the constants. This gives the following theorem for general numerical ranges.

5 Weyl calculus for time-ordered evolution↩︎

The two-dimensional Weyl calculus developed in 3 applies directly to the function \(f(z)=e^{-z}=e^{-(x+\mathrm{i}y)}\) on \(\Omega=\{z:\mathrm{Re}(z)\ge0\}\). In 6, we may choose \(\hat{h}=\sqrt{2\pi}\delta_1(v)r(u)\,\mathrm{d}u\) as a product of two one-dimensional measures. Here \(\delta_1(v)\) is the Dirac measure supported at \(1\), and \(r(u)\,\mathrm{d}u\) is an absolutely continuous measure with inverse Fourier transform \(\check{r}(x)\approx e^{-x}\) for \(x\ge0\). However, 6 does not directly apply to the time-dependent ODE \[\label{eq:32time32dependent32ODE} \frac{\,\mathrm{d}}{\,\mathrm{d}t} \vec{\psi}(t) = -A(t)\vec{\psi}(t)\tag{46}\] which is an important application of LCHS. In this section, we show how to adapt the Weyl-calculus framework to the time-ordered operator \[\label{eq:32time32ordered32operator} \mathcal{T}e^{-\int_0^t A(s) \,\mathrm{d}s}\tag{47}\] that solves 46 .

Let \(L(s)\) and \(H(s)\) be Hermitian matrix-valued functions. Suppose a function \(f\) admits the Fourier–Stieltjes representation \[\label{eq:321d32FS32f} f(x) = (2 \pi)^{-1/2} \int_{\mathbb{R}} e^{-\mathrm{i}ux} \,\mathrm{d}\mu(u)\tag{48}\] for some finite complex Borel measure \(\mu\). We define \[\label{eq:32def32whfl} \mathrm{Op}_{H}(f)\left(L\right):=(2 \pi)^{-1/2} \int_{\mathbb{R}} \mathcal{T}e^{-\mathrm{i}\int_0^t (u L(s)+H(s))\,\mathrm{d}s} \,\mathrm{d}\mu(u).\tag{49}\]

5.1 Improved prefactor for dissipative dynamics simulation↩︎

In this subsection, we present a new kernel function for the LCHS simulation of \[\mathcal{T}e^{-\int_0^t A(s) \,\mathrm{d}s}\] for a time-dependent matrix \(A(t)\). This is the first and most extensively studied application of LCHS. Ref. [9] gives an efficient kernel construction using a contour-integral argument; here we construct a kernel with a better prefactor.

For positive parameters \(R\) and \(C\), define \[\begin{align} \hat{g}_{R,C}(u) = \sqrt{\frac{2}{\pi}} e^{2C \arctan\left(\frac{1}{R}\right)} e^{\mathrm{i}C \log\left(\frac{R-u}{R+u}\right)} \frac{1}{1+u^2} \end{align}\] on the interval \((-R,R)\), and set \[\begin{align} \label{eq:32fRC32definition} f_{R,C}(x) := \frac{1}{\sqrt{2\pi}}\int_{-R}^R\hat{g}_{R,C}(u) e^{-\mathrm{i}ux} \,\mathrm{d}u. \end{align}\tag{53}\]

Combining the preceding lemma with 10 and optimizing the parameters \(R\) and \(C\) gives the following theorem.

For comparison, the state-of-the-art scaling from Ref. [9] is asymptotically \(2.57\log_2(\frac{1}{\epsilon}) \approx 3.71\ln\frac{1}{\epsilon}\). Our explicit parameter choice in the new kernel yields \(\frac{2e}{\pi}\ln\frac{1}{\epsilon} \approx 1.73\ln\frac{1}{\epsilon}\), giving a \(2.1\times\) improvement in the asymptotic constant.

An even tighter cost estimate is obtained by numerically optimizing 57 ; the results are shown in 1.

6 Cost optimization of discrete LCHS formulas↩︎

The continuous LCHS formula ?? is not directly implementable in a coherent quantum circuit 7.1. Previous works [4]–[6], [9], [16] discretize such formulas using quadrature rules, and the same approach applies to the explicit continuous constructions in this paper. However, quadrature-based discretization can be inefficient in two ways:

1. If the weight function \(\hat{h}\) is not sufficiently smooth, or if the quadrature rule is not spectrally convergent, the number of discretization terms can be \(\mathcal{O}\left(\mathop{\mathrm{poly}}(\epsilon^{-1})\right)\). In the LCU implementation, this leads to an additional \(\mathcal{O}\left(\mathop{\mathrm{polylog}}(\epsilon^{-1})\right)\) ancilla-qubit cost in the \(\mathrm{PREP}\) operators 70 .

2. The cost \(\Lambda\) associated with an explicit construction \(\hat{h}(u,v)\) may have optimal asymptotic scaling but suboptimal constants. A direct numerical optimization method can give more favorable resource estimates.

For the exponential function \(f(z)=e^{-z}\), the works [6], [9] address the second issue by optimizing over parametrized families of weight functions. The output, however, is still a continuous kernel and must be discretized afterward. Moreover, the parametrizations are chosen from intuition, which may not be optimal and may be difficult to generalize beyond \(e^{-z}\).

The Weyl-calculus framework developed in this paper suggests a direct alternative: optimize the weights of the discrete LCHS formula themselves, without first choosing a continuous ansatz.

Assume that the truncation radius \(R\) and total weight budget \(\lambda\) are fixed, and choose a grid \[\mathcal{G} = \{(u_j, v_j)\}_{j=1,\ldots,J}\subset D_R.\] We seek weights \(\mathbf{w}=(w_j)_{1\le j\le J}\) solving the convex optimization problem \[\label{eq:32general32opt32problem} \begin{align} \epsilon = \min_{\mathbf{w}} &\; \mathchoice {\left\lVert f(x+\mathrm{i}y) - \sum_{j=1}^J w_j\, e^{-\mathrm{i}(u_jx+v_jy)}\right\rVert} {\lVert f(x+\mathrm{i}y) - \sum_{j=1}^J w_j\, e^{-\mathrm{i}(u_jx+v_jy)}\rVert} {\lVert f(x+\mathrm{i}y) - \sum_{j=1}^J w_j\, e^{-\mathrm{i}(u_jx+v_jy)}\rVert} {\lVert f(x+\mathrm{i}y) - \sum_{j=1}^J w_j\, e^{-\mathrm{i}(u_jx+v_jy)}\rVert} _{W(\Omega)}\\ &\text{s.t.}\quad \sum_{j=1}^J|w_j|\le \lambda \end{align}\tag{58}\] where \(\Omega\) is a region containing the numerical range of \(X+\mathrm{i}Y\). The function \(h(x,y)=\sum_{j=1}^J w_j e^{-\mathrm{i}(u_jx+v_jy)}\) is the inverse Fourier transform of the measure \(\mu=2\pi\sum_{j=1}^J w_j\delta_{(u_j,v_j)}\). By 6, \[\mathchoice {\left\lVert f(X+\mathrm{i}Y) - \sum_{j=1}^J w_j\, e^{-\mathrm{i}(u_jX+v_jY)}\right\rVert} {\lVert f(X+\mathrm{i}Y) - \sum_{j=1}^J w_j\, e^{-\mathrm{i}(u_jX+v_jY)}\rVert} {\lVert f(X+\mathrm{i}Y) - \sum_{j=1}^J w_j\, e^{-\mathrm{i}(u_jX+v_jY)}\rVert} {\lVert f(X+\mathrm{i}Y) - \sum_{j=1}^J w_j\, e^{-\mathrm{i}(u_jX+v_jY)}\rVert} \le \epsilon,\] and the cost is \(\Lambda=\lambda R\). In practice, one may instead fix \(\Lambda\) and optimize the resulting error \(\epsilon\). This can be done by trying different values of \(R\), setting \(\lambda=\Lambda/R\), and choosing the smallest certified error.

Computing the exact local Wiener norm is itself difficult. Therefore, when implementing 58 , we replace it by explicit upper bounds. For common domains such as intervals, squares, and discs, the resulting convex programs are described below.

6.1 Convex optimization on \([0,T]\): application to dissipative dynamics↩︎

To obtain a discrete LCHS formula for \(\mathcal{T}e^{-\int_0^t A(s)\,\mathrm{d}s}\) under the dissipativity condition \(L(s)=\frac{A(s)+A(s)^\dagger}{2}\succeq0\), 9 suggests choosing \(T\ge \mathchoice {\left\lVert L\right\rVert} {\lVert L\rVert} {\lVert L\rVert} {\lVert L\rVert} _{L^1([0,t])}=\int_0^t \mathchoice {\left\lVert L(s)\right\rVert} {\lVert L(s)\rVert} {\lVert L(s)\rVert} {\lVert L(s)\rVert} \,\mathrm{d}s\) and approximating \(e^{-x}\) on \([0,T]\). The same lemma also applies to non-dissipative dynamics by approximating \(e^{-x}\) on the interval \([\int_0^t l(s)\,\mathrm{d}s,\int_0^t r(s)\,\mathrm{d}s]\); this extension is straightforward, so we omit it for conciseness.

Let \(s=\frac{\pi}{T}\), and choose a truncation interval \([-R,R]\) such that \(N=R/s\) is an integer. We use the finite Fourier expansion \[S(x)=s\sum_{j=-N}^{N}g_j e^{-\mathrm{i}u_jx}, \qquad u_j=js,\] to approximate \(e^{-x}\) for \(x\in[0,T]\), where the coefficients \(g_j\in\mathbb{C}\) are variables to be determined. The objective is to minimize the local Wiener norm of the residual \[\eta(x)=S(x)-e^{-x},\qquad x\in[0,T].\] Denote the cosine expansion of \(\eta\) on \([0,T]\) by \[\eta(x)=\frac{a_0}{2}+\sum_{m=1}^{\infty}a_m\cos(c_mx), \qquad c_m=\frac{m\pi}{T}.\] Then \[\mathchoice {\left\lVert\eta\right\rVert} {\lVert\eta\rVert} {\lVert\eta\rVert} {\lVert\eta\rVert} _{W([0,T])}\le \frac{|a_0|}{2}+\sum_{m=1}^{\infty}|a_m|.\] We use the cosine expansion coefficients to upper bound the local Wiener norm because they are absolutely summable for a smooth function on a finite interval. In particular, for a fixed smooth residual \(\eta\), the coefficients \(a_m\) decay at least quadratically in \(m\).

The coefficients \(a_m\) can be computed analytically. Define \[\begin{align} D_{j,m} &=\frac{2}{T}\int_0^T e^{-\mathrm{i}u_jx}\cos(c_mx)\,\mathrm{d}x \notag\\ &=\frac{1}{T}\left( \frac{1-e^{-\mathrm{i}(u_j-c_m)T}}{\mathrm{i}(u_j-c_m)} +\frac{1-e^{-\mathrm{i}(u_j+c_m)T}}{\mathrm{i}(u_j+c_m)} \right),\label{eq:32Djm32exp32showcase} \end{align}\tag{59}\] where the values at removable singularities are clear. Also set \[C_m=\frac{2}{T}\int_0^T e^{-x}\cos(c_mx)\,\mathrm{d}x =\frac{2/T}{1+c_m^2}\left[1-(-1)^m e^{-T}\right].\] Thus \[a_m=s\sum_{j=-N}^{N}g_jD_{j,m}-C_m .\]

For fixed \(R\) and \(\lambda\), truncating the cosine series at a sufficiently large value \(M\) gives the finite convex program \[\label{eq:32exp32cvx32constraint} \epsilon_M= \min_{\{g_j\}_{j=-N}^{N}}\; \frac{1}{2}\left|s\sum_{j=-N}^{N}g_jD_{j,0}-C_0\right| +\sum_{m=1}^{M}\left|s\sum_{j=-N}^{N}g_jD_{j,m}-C_m\right|\tag{60}\] subject to \[s\sum_{j=-N}^{N}|g_j|\le \lambda .\] This is a second-order cone program (SOCP), since both the objective and the constraint are sums of moduli of complex affine functions. It can be solved efficiently using standard convex optimization solvers. To verify a solution, we use a larger value \(M_{\text{test}}\) in the right-hand side of 60 , giving a sharper posterior error bound.

The practical implementation imposes the conjugate symmetry \(\overline{g_{-j}}=g_j\), since the target function \(e^{-x}\) is real-valued. This reduces the number of variables by roughly a factor of two.

We run this solver using the \(R\) and \(\lambda\) values from 1. The resulting \(\epsilon\) values, shown in 2 (a), are almost identical to those obtained from the continuous distribution \(\hat{g}_{R,C}\), confirming the effectiveness of the discrete optimization. These results use \(T=50\), and we also verify that they are insensitive to this choice. We use \(M_{\text{opt}}=3000\) in the optimization and certify the residual error with 60 using \(M_{\text{test}}=20000\).

2 (b) shows the optimized discrete weights \(g_j\) for \(T=10\), \(R=6.90\), and \(\lambda=2.28\). The weights are close to the values \(\frac{1}{\sqrt{2\pi}}\hat{g}_{R,C}(u_j)\). In fact, the closed-form kernel \(\hat{g}_{R,C}\) was originally inspired by this optimization result. Even with the explicit formula for \(\hat{g}_{R,C}\), the optimized discrete weights remain useful for implementation: they live on a grid of size only \(\mathcal{O}\left(R\right)=\mathcal{O}\left(\log\frac{1}{\epsilon}\right)\). By comparison, discretizing the continuous integral 53 by quadrature may require \(\mathcal{O}\left(\mathop{\mathrm{poly}}(\frac{1}{\epsilon})\right)\) terms because \(\hat{g}_{R,C}\) has singularities near \(\pm R\).

6.2 Convex optimization on \([-1,1]^2\): application to matrix square roots↩︎

If the numerical range of \(A\) is contained in the square \(\Omega = [-1,1]^2\), we can approximate \(f(x+\mathrm{i}y)\) using \[S(x,y) = \sum_{(j,k)\in\mathcal{G}_R} w_{j,k} e^{-\mathrm{i}\frac{\pi}{2}(jx+ky)},\] on \(\Omega\), where \[\label{eq:32disc32frequency32set} \mathcal{G}_R =\left\{(j,k)\in\mathbb{Z}^2: \left(\frac{\pi}{2}\right)^2(j^2+k^2)\le R^2\right\}.\tag{61}\] The optimization problem is similar to the one-dimensional case, and the local Wiener norm can again be upper-bounded using a cosine expansion. Suppose that the residual \(\eta(x,y)=S(x,y)-f(x+\mathrm{i}y)\) has the cosine expansion \[\eta(x,y) = \sum_{m,n=0}^\infty a_{m,n} \cos\left(\frac{m\pi}{2}(x+1)\right) \cos\left(\frac{n\pi}{2}(y+1)\right).\] Then \[\mathchoice {\left\lVert\eta\right\rVert} {\lVert\eta\rVert} {\lVert\eta\rVert} {\lVert\eta\rVert} _{W(\Omega)} \le \sum_{m,n=0}^\infty |a_{m,n}|.\] The coefficients \(a_{m,n}\) can be computed analytically by taking inner products with the basis functions. In particular, set \[\begin{align} D_{m,n}^{j,k} &= \frac{1}{(1+\delta_{m,0})(1+\delta_{n,0})}\int_{-1}^1\int_{-1}^1 e^{-\mathrm{i}\frac{\pi}{2}(jx+ky)} \cos\left(\frac{m\pi}{2}(x+1)\right) \cos\left(\frac{n\pi}{2}(y+1)\right) \mathrm{d}x \mathrm{d}y, \\ C_{m,n} &= \frac{1}{(1+\delta_{m,0})(1+\delta_{n,0})}\int_{-1}^1\int_{-1}^1 f(x+\mathrm{i}y) \cos\left(\frac{m\pi}{2}(x+1)\right) \cos\left(\frac{n\pi}{2}(y+1)\right) \mathrm{d}x \mathrm{d}y. \end{align}\] The resulting optimization problem is \[\begin{align} \epsilon_{M} = \min_{\{w_{j,k}\}} &\;\sum_{m=0}^M \sum_{n=0}^M \left| \sum_{(j,k)\in\mathcal{G}_R} w_{j,k} D_{m,n}^{j,k} - C_{m,n} \right| \\ &\text{s.t.}\quad \sum_{(j,k)\in\mathcal{G}_R} |w_{j,k}| \le \lambda, \end{align}\] As in the one-dimensional case, we use a sufficiently large \(M_{\text{opt}}\) during optimization and a larger \(M_{\text{test}}\) for posterior verification.

We apply this optimization approach to compute LCHS weights for the matrix square root \(\sqrt{A}\), assuming that the numerical range of \(A\) is contained in the open right half-plane \(\{\mathrm{Re}(z)>0\}\). To adapt the problem to the domain \(\Omega=[-1,1]^2\), we apply a linear transformation \(A'=\beta A-(1+\rho)I\) so that the numerical range of \(A'\) is contained in \(\Omega\). Thus we focus on functions \[f(z) = \sqrt{(1+\rho)+z}\] with \(\rho>0\), and approximate them on \(\Omega=[-1,1]^2\).

We present numerical results for \(\rho=1,0.5,0.25,0.125\). For a fixed budget \(\Lambda\), we try different values of \(R\), set \(\lambda=\Lambda/R\), run the optimization with \(M_{\text{opt}}=60\) and \(M_{\text{test}}=1000\), and report the smallest certified \(\epsilon\). The results are shown in 3.

It is useful to compare these numerical results with the theoretical bounds. The bound in 7 applies because \(D_1\subset[-1,1]^2\), so controlling the square residual also controls the disc residual. The parameter \(\rho\) has the same meaning here because \(\sqrt{(1+\rho)+z}\) has finite local Wiener norm on \(D_{1+\rho}\). In fact, \(\sqrt{(1+\rho)+z}\) can be extended to an \(H^{1.25}(\mathbb{R}^2)\) function, and one can prove \(\mathchoice {\left\lVert g\right\rVert} {\lVert g\rVert} {\lVert g\rVert} {\lVert g\rVert} _{W(\mathbb{R}^2)}\le C \mathchoice {\left\lVert g\right\rVert} {\lVert g\rVert} {\lVert g\rVert} {\lVert g\rVert} _{H^{1.25}(\mathbb{R}^2)}\) for some constant \(C\) by replacing the \((1+u^2+v^2)^2\) factor in 9 with \((1+u^2+v^2)^{1.25}\). A simple lower bound is \(\mathchoice {\left\lVert\sqrt{(1+\rho)+z}\right\rVert} {\lVert\sqrt{(1+\rho)+z}\rVert} {\lVert\sqrt{(1+\rho)+z}\rVert} {\lVert\sqrt{(1+\rho)+z}\rVert} _{W(D_{1+\rho})}\ge \mathchoice {\left\lVert\sqrt{(1+\rho)+z}\right\rVert} {\lVert\sqrt{(1+\rho)+z}\rVert} {\lVert\sqrt{(1+\rho)+z}\rVert} {\lVert\sqrt{(1+\rho)+z}\rVert} _{L^{\infty}(D_{1+\rho})}\ge\sqrt{2}\). Therefore, the results in 3 have lower cost than even the most optimistic upper bound obtainable from 7, namely \(\Lambda \sim \frac{3\sqrt{2}}{\rho}\ln \frac{30\sqrt{2}}{\rho\epsilon}\).

6.3 Convex optimization on \(D_1\): example of \(f(z)=z^4\)↩︎

If a function has moderate magnitude on the unit disc \(D_1\) but large values on the rest of the square \([-1,1]^2\), it is wasteful to measure the residual by a tensor-product cosine expansion on the whole square. In this case, we instead optimize the residual directly on \(D_1\). We still approximate \(f(x+\mathrm{i}y)\) on \(D_1\) by \[\label{eq:32disc32trig32approximation} S(x,y) = \sum_{(j,k)\in\mathcal{G}_R} w_{j,k} e^{-\mathrm{i}\frac{\pi}{2}(jx+ky)},\tag{62}\] where \(w_{j,k}\in\mathbb{C}\).

We use the Neumann Fourier–Bessel basis on the disc to measure the residual. Write \(x=r\cos\theta\) and \(y=r\sin\theta\). Let \(J_p\) be the Bessel function of the first kind, and let \(\beta_{p,\ell}\) be the \(\ell\)-th non-negative zero of \(J_p'\). Thus \(\beta_{0,1}=0\), while otherwise \(\beta_{p,\ell}>0\). This gives the following unnormalized orthogonal basis for \(L^2(D_1)\): \[\label{eq:32disc32FB32basis} \Psi_{m,\ell}(r,\theta) =J_{|m|}(\beta_{|m|,\ell}r)e^{\mathrm{i}m\theta}, \qquad m\in\mathbb{Z},\quad \ell=1,2,\ldots .\tag{63}\] The mode \((m,\ell)=(0,1)\) is the constant function. The zero convention, the Neumann boundary condition, and the orthogonality relation are recalled in 10.1. The normalization constant is \[\label{eq:32disc32FB32norm32constant} (2\pi)N_{p,\ell}= \mathchoice {\left\lVert\Psi_{m,\ell}\right\rVert} {\lVert\Psi_{m,\ell}\rVert} {\lVert\Psi_{m,\ell}\rVert} {\lVert\Psi_{m,\ell}\rVert} _{L^2(D_1)}^2 = (2\pi)\int_0^1J_p(\beta_{p,\ell}r)^2r\,\mathrm{d}r.\tag{64}\] For the residual \[\eta(r,\theta)=S(r\cos\theta,r\sin\theta)-f(re^{\mathrm{i}\theta}),\] we expand it as \[\label{eq:32bessel32a95ml} \eta=\sum_{m\in\mathbb{Z}}\sum_{\ell\ge1}a_{m,\ell}\Psi_{m,\ell}.\tag{65}\] In 10.2, we establish \(\mathchoice {\left\lVert\Psi_{m,\ell}\right\rVert} {\lVert\Psi_{m,\ell}\rVert} {\lVert\Psi_{m,\ell}\rVert} {\lVert\Psi_{m,\ell}\rVert} _{W(D_1)}\le 1\). Consequently, \[\mathchoice {\left\lVert\eta\right\rVert} {\lVert\eta\rVert} {\lVert\eta\rVert} {\lVert\eta\rVert} _{W(D_1)} \le \sum_{m\in\mathbb{Z}}\sum_{\ell\ge1}|a_{m,\ell}|.\] This is the disc analogue of the cosine coefficient bound used in the one-dimensional and square cases. It remains to write the residual coefficients as affine functions of the unknown weights. The Neumann Fourier–Bessel expansion of the plane wave \(e^{-\mathrm{i}\frac{\pi}{2}(jx+ky)}\) is \[e^{-\mathrm{i}\frac{\pi}{2}(jx+ky)} = \sum_{m\in\mathbb{Z}}\sum_{\ell\ge1}P_{m,\ell}^{j,k}\Psi_{m,\ell}(r,\theta),\] where the coefficient \(P_{m,\ell}^{j,k}\) has a closed-form expression recorded in 10.3. The target coefficients are \[\label{eq:32disc32target32coefficient} C_{m,\ell} =\frac{1}{2\pi N_{|m|,\ell}} \int_0^{2\pi}\int_0^1 f(re^{\mathrm{i}\theta})J_{|m|}(\beta_{|m|,\ell}r)e^{-\mathrm{i}m\theta}\,r\,\mathrm{d}r\,\mathrm{d}\theta.\tag{66}\] Therefore, the residual coefficient is the affine expression \[\label{eq:32disc32residual32coefficient32affine} a_{m,\ell} =\sum_{(j,k)\in\mathcal{G}_R}w_{j,k}P_{m,\ell}^{j,k}-C_{m,\ell}.\tag{67}\]

For computation, choose an integer \(M\) and truncate to \[\label{eq:32disc32FB32truncation32set} \mathcal{I}_M =\left\{(m,\ell): \beta_{|m|,\ell}\le \frac{\pi}{2}M\right\}.\tag{68}\] For frequency radius \(R\) and coefficient budget \(\lambda\), solve \[\label{eq:32disc32FB32cvx32program} \begin{align} \epsilon_M= \min_{\{w_{j,k}\}_{(j,k)\in\mathcal{G}_R}} &\quad \sum_{(m,\ell)\in\mathcal{I}_M} \left|\sum_{(j,k)\in\mathcal{G}_R}w_{j,k}P_{m,\ell}^{j,k}-C_{m,\ell}\right| \\ \text{s.t.} &\quad \sum_{(j,k)\in\mathcal{G}_R}|w_{j,k}|\le \lambda . \end{align}\tag{69}\] This is again a second-order cone program. After obtaining the optimizer, one can recompute the coefficient sum on a larger set \(\mathcal{I}_{M_{\mathrm{test}}}\) and use the enlarged sum as the reported disc error. This verifier, however, is time-consuming even for moderately large \(M_{\mathrm{test}}\), because we do not have a simple FFT analogue for computing the Neumann Fourier–Bessel coefficients. In practice, we use a more efficient verifier, discussed in 10.4. The idea is similar to the derivative bounds 9 and 10 , but the practical verifier gives a much tighter upper bound through more refined calculations.

For functions that become large on \([-1,1]^2\setminus D_1\), the disc optimizer is preferable to the square optimizer in 6.2 whenever we know \(\mathchoice {\left\lVert A\right\rVert} {\lVert A\rVert} {\lVert A\rVert} {\lVert A\rVert} \le1\). We illustrate this with \(f(z)=z^4\). This example is only illustrative, since matrix monomials can be computed by more direct quantum methods. We fix \(\lambda=12\) and \(M_{\text{opt}}=75\) for both the square and disc optimizers, and compare the truncation radius \(R\) required for different values of \(\epsilon\) in 4. The disc optimizer has a significantly lower cost, as expected.

We also remark that the Neumann Fourier–Bessel coefficient \(\ell^1\) norm is usually not an optimal estimate of the local Wiener norm. Therefore, for functions that remain moderate on \([-1,1]^2\setminus D_1\), the square optimizer may give comparable or even better results.

7 Quantum implementation of LCHS↩︎

Given an LCHS formula of the form ?? or 4 , we briefly discuss how to implement it on a quantum computer. The goal of this section is to explain how the quantum cost depends on the LCHS parameters. The material is standard in LCHS-based methods; see [4], [5], [9], [16] for more details.

There are two common ways of implementing LCHS formulas. The first is a coherent implementation, which combines LCU [33] with Hamiltonian simulation to block-encode the entire linear combination. The Hamiltonian simulation subroutine can be implemented using QSVT [3], Trotter formulas [34], or randomized methods [35], [36]. This is the natural choice when the goal is to prepare a state proportional to \(f(A)\ket{\psi}\) or to use \(f(A)\) as a subroutine in a larger quantum algorithm. The second is a hybrid implementation, designed for estimating observables such as \(\bra{\psi} f(A)^\dagger O f(A)\ket{\psi}\). It samples terms from the LCHS formula classically and evaluates the corresponding quantum correlation functions one at a time. The coherent approach is usually better in asymptotic quantum resource scaling, while the hybrid approach is more suitable for near-term devices: it uses fewer ancilla qubits, avoids a large coherently controlled select oracle, and yields shallower circuits.

7.1 Coherent implementation↩︎

The coherent implementation is most directly stated for the discrete formula 4 , which we discuss in this subsection. We begin with the standard block-encoding model for representing matrices on quantum computers.

The LCU method [33] is a standard technique for implementing matrix sums. Given unitaries \(\{U_j\}\) and coefficients \(\{w_j\}\), LCU implements a \(\lambda\)-block-encoding of \(\sum_j w_jU_j\), where \(\lambda=\sum_j |w_j|\). In detail, one uses the prepare operators \[\label{eq:32prep32operators} \mathrm{PREP}_l\ket{0} = \frac{1}{\sqrt{\lambda}}\sum_{j=1}^J \sqrt{|w_j|}\ket{j},\qquad \mathrm{PREP}_r\ket{0} = \frac{1}{\sqrt{\lambda}}\sum_{j=1}^J \sqrt{|w_j|}e^{\mathrm{i}\mathop{\mathrm{Arg}}(w_j)}\ket{j},\tag{70}\] and the select operator \[\mathrm{SEL} = \sum_{j=1}^J \ket{j}\bra{j}\otimes U_j.\] Then the product \[(\mathrm{PREP}_l^{\dagger}\otimes I)\;\mathrm{SEL}\;(\mathrm{PREP}_r\otimes I)\] is a \(\lambda\)-block-encoding of \(\sum_j w_jU_j\).

Suppose \(U_A\) and \(U_{A^\dagger}\) are \(\alpha\)-block-encodings of \(A\) and \(A^\dagger\), respectively, and let \(u,v\in\mathbb{R}\). Applying LCU to \(U_A\) and \(U_{A^\dagger}\) gives an \((\alpha\sqrt{u^2+v^2})\)-block-encoding of \(uX+vY\), since \[uX+vY = \frac{u-\mathrm{i}v}{2}A + \frac{u+\mathrm{i}v}{2}A^\dagger,\] and \(\frac{|u-\mathrm{i}v|}{2}+\frac{|u+\mathrm{i}v|}{2}=\sqrt{u^2+v^2}\). We can then apply QSVT-based Hamiltonian simulation [3], [8] to implement \(e^{-\mathrm{i}(uX+vY)}\) with query complexity \(\mathcal{O}\left(\alpha\sqrt{u^2+v^2}+\log(\epsilon^{-1})\right)\).

Finally, applying LCU again to the discrete LCHS formula 4 , with \(U_j=e^{-\mathrm{i}(u_jX+v_jY)}\), gives a \(\lambda\)-block-encoding of \(f(A)\) with \(\lambda=\sum_j |w_j|\). The query complexity of implementing the select operator \[\mathrm{SEL} = \sum_{j=1}^J \ket{j}\bra{j}\otimes e^{-\mathrm{i}(u_jX+v_jY)}\] is the maximum query complexity of implementing \(e^{-\mathrm{i}(u_jX+v_jY)}\) over all \(j\), namely \(\mathcal{O}\left(\alpha R+\log(\epsilon^{-1})\right)\). This is because the select operator can be implemented by a single Hamiltonian simulation circuit with the index register as a control. This procedure can be illustrated as \[\begin{align} \text{controlled access to } U_A,U_{A^\dagger} &\xrightarrow{\text{controlled version of LCU}} \sum_{j=1}^J \ket{j}\bra{j}\otimes \begin{bmatrix} \frac{u_jX+v_jY}{\alpha R} & \cdot \\ \cdot & \cdot \end{bmatrix}\\ &\xrightarrow{\text{QSVT Hamiltonian simulation}} \sum_{j=1}^J \ket{j}\bra{j}\otimes e^{-\mathrm{i}(u_jX+v_jY)} = \mathrm{SEL}, \end{align}\] and we refer readers to [4], [9] for circuit details.

If we want to amplify the block-encoding normalization from \(\lambda\) toward \(\mathchoice {\left\lVert f(A)\right\rVert} {\lVert f(A)\rVert} {\lVert f(A)\rVert} {\lVert f(A)\rVert}\), we can use uniform singular value amplification [3]; the query complexity then increases by a factor of \(\lambda/ \mathchoice {\left\lVert f(A)\right\rVert} {\lVert f(A)\rVert} {\lVert f(A)\rVert} {\lVert f(A)\rVert}\). Thus the final query complexity is proportional to \(\frac{\lambda}{ \mathchoice {\left\lVert f(A)\right\rVert} {\lVert f(A)\rVert} {\lVert f(A)\rVert} {\lVert f(A)\rVert} }(\alpha R+\log(\epsilon^{-1}))\). Ignoring constants and using \(R=\Omega(\log(\epsilon^{-1}))\), the query complexity in \(A\) and \(A^\dagger\) is proportional to \[\Lambda := R \lambda.\]

If the goal is to prepare the normalized state proportional to \(f(A)\ket{\psi}\) for a given \(\ket{\psi}\), amplitude amplification [10], [37] also gives the same proportional dependence on \(\Lambda\). This justifies using \(\Lambda\) as the cost function of an LCHS formula.

7.2 Hybrid implementation↩︎

The hybrid implementation applies when the goal is only to evaluate an observable. Suppose \(O\) is Hermitian and the desired quantity is \(\bra{\psi} f(A)^\dagger O f(A)\ket{\psi}\). Using the discrete approximation \(f(A)=\sum_j w_jU_j\), this observable becomes \[\label{eq:32hybrid32observable32expansion} \bra{\psi} f(A)^\dagger O f(A)\ket{\psi} = \sum_{j,k=1}^J \overline{w_j}w_k \bra{\psi}U_j^\dagger O U_k\ket{\psi}.\tag{71}\] The hybrid method samples an index pair \((j,k)\) with probability proportional to \(|w_jw_k|\), and then estimates the matrix element \(o_{j,k}:=\bra{\psi}U_j^\dagger O U_k\ket{\psi}\) using the Hadamard test for non-unitary matrices [19] and amplitude estimation [10]. Define \(Z=e^{\mathrm{i}\mathop{\mathrm{Arg}}(\overline{w_j}w_k)}o_{j,k}\). With sufficiently many samples, the sample average converges to \(\mathbb{E}Z=\lambda^{-2}\bra{\psi}f(A)^\dagger O f(A)\ket{\psi}\), which is the desired observable up to the normalization factor \(\lambda^2\).

A direct bounded-variance estimate [4] gives sample complexity proportional to \(\lambda^4 \mathchoice {\left\lVert O\right\rVert} {\lVert O\rVert} {\lVert O\rVert} {\lVert O\rVert} ^2/\epsilon^2\), up to logarithmic factors in the failure probability, for additive error \(\epsilon\). Since each sample requires Hamiltonian simulation for time \(R\), the total query complexity is proportional to \(R\lambda^4 \mathchoice {\left\lVert O\right\rVert} {\lVert O\rVert} {\lVert O\rVert} {\lVert O\rVert} ^2/\epsilon^2\). In this setting, however, the circuit depth \(R\) is often the most important resource, because hybrid implementations are usually considered for early fault-tolerant devices with limited coherence time. The same procedure also applies directly to the continuous formula ?? , since one can sample \((u,v)\) from a continuous distribution.

8 Conclusion and discussion↩︎

This work shows that Weyl calculus provides a natural language for LCHS beyond the matrix exponential. The main technical advantage is that the noncommutativity of the Hermitian and skew-Hermitian parts of a matrix is absorbed into the Weyl quantization map, while the approximation problem remains a scalar approximation problem on the numerical range. In this sense, 6 separates two issues that were intertwined in earlier LCHS constructions: finding a good scalar approximation to \(f(x+\mathrm{i}y)\) and implementing the corresponding plane-wave expansion by Hamiltonian simulation. The framework also treats discrete formulas as first-class objects through the Fourier–Stieltjes viewpoint, which is important for coherent quantum implementation because the final algorithm ultimately uses a finite LCU expansion.

A first direction for future work is to make the numerical optimization framework more efficient. The convex program 58 directly optimizes the discrete weights and therefore avoids both quadrature error and the choice of a parametrized ansatz. Its practical performance, however, depends on how tightly and efficiently one can upper bound the local Wiener norm of the residual. The cosine and Fourier–Bessel coefficient bounds used in this paper are convenient and convex, but they are not expected to be optimal. It would be useful to develop sharper verifiers for common domains and faster structured solvers for the resulting SOCPs.

For non-unitary dissipative dynamics simulation, our kernel improves the asymptotic prefactor of the LCHS cost to \(\frac{2e}{\pi}\ln(1/\epsilon)\), which is within about a factor of \(1.5\) of the lower bound from [9]. Closing this gap and matching the upper and lower bounds remains an interesting open problem.

Another question is whether Weyl-calculus-based LCHS can lead to quantum algorithms beyond eigenvalue transformation. The time-dependent dissipative linear ODE example suggests that this may be possible. It would be interesting to explore whether the framework can be adapted to other differential equations or linear algebra tasks.

9 Wiener norm and Fourier–Stieltjes algebra↩︎

9.1 The Fourier–Stieltjes algebra \(B(\mathbb{R}^n)\)↩︎

Let \(M(\mathbb{R}^n)\) be the Banach space of finite complex Borel measures on \(\mathbb{R}^n\), equipped with the total variation norm. The Fourier–Stieltjes algebra \(B(\mathbb{R}^n)\) consists of all functions of the form \[\label{eq:32B32representation32appendix} g(\mathbf{x}) = (2 \pi)^{-n / 2} \int_{\mathbb{R}^n} e^{-\mathrm{i}\boldsymbol{\xi}\cdot \mathbf{x}} \,\mathrm{d}\mu(\boldsymbol{\xi}), \qquad \mu\in M(\mathbb{R}^n).\tag{72}\] The Fourier transform of a finite Borel measure is unique. Thus, if \(\hat{g}=\mu\), we define \[\label{eq:32B32norm32appendix} \mathchoice {\left\lVert g\right\rVert} {\lVert g\rVert} {\lVert g\rVert} {\lVert g\rVert} _W := (2 \pi)^{-n / 2} |\mu|(\mathbb{R}^n).\tag{73}\] In addition to the absolutely continuous case \(\mu=\hat{g}(\boldsymbol{\xi})\,\mathrm{d}\boldsymbol{\xi}\), an important example is the Dirac comb \[\mu = \sum_{\mathbf{k}\in\mathbb{Z}^n} c_{\mathbf{k}} \delta_{a\mathbf{k}},\] whose inverse Fourier transform is the Fourier series \[g(\mathbf{x}) = (2\pi)^{-n/2}\sum_{\mathbf{k}\in\mathbb{Z}^n} c_{\mathbf{k}} e^{-\mathrm{i}a\mathbf{k}\cdot \mathbf{x}}.\]

9.2 Local Wiener norm↩︎

We now prove 8 , which states that the local Wiener norm on a compact set \(\Omega\) can also be computed using extensions in \(A(\mathbb{R}^n)\) whose Fourier transforms are smooth.

Let \(\Omega\) be a compact subset of \(\mathbb{R}^n\). For any \(\epsilon > 0\), [38] gives a function \(\chi\in A(\mathbb{R}^n)\) such that \(\chi=1\) on \(\Omega\), \(\chi\) has compact support, and \(\mathchoice {\left\lVert\chi\right\rVert} {\lVert\chi\rVert} {\lVert\chi\rVert} {\lVert\chi\rVert} _W=(2\pi)^{-n/2} \mathchoice {\left\lVert\hat{\chi}\right\rVert} {\lVert\hat{\chi}\rVert} {\lVert\hat{\chi}\rVert} {\lVert\hat{\chi}\rVert} _{L^1(\mathbb{R}^n)}\le 1+\epsilon\). Since \(\chi\) has compact support, \(\hat{\chi}\) is smooth. For any extension \(\tilde{g}\in B(\mathbb{R}^n)\) of \(g\), its Fourier transform \((\tilde{g})^\wedge = \tilde{\mu}\) is a finite measure. Consider the product \(h = \chi \tilde{g}\). Under our Fourier convention, \[\hat{h} = (2\pi)^{-n/2}\hat{\chi}*\tilde{\mu}\in C^{\infty}(\mathbb{R}^n),\] because \(\hat{\chi}\) is smooth and \(\tilde{\mu}\) is finite. Moreover, Young’s inequality gives \[\mathchoice {\left\lVert\hat{h}\right\rVert} {\lVert\hat{h}\rVert} {\lVert\hat{h}\rVert} {\lVert\hat{h}\rVert} _{L^1(\mathbb{R}^n)} \le (2\pi)^{-n/2} \mathchoice {\left\lVert\hat{\chi}\right\rVert} {\lVert\hat{\chi}\rVert} {\lVert\hat{\chi}\rVert} {\lVert\hat{\chi}\rVert} _{L^1(\mathbb{R}^n)} \mathchoice {\left|\tilde{\mu}\right|} {|\tilde{\mu}|} {|\tilde{\mu}|} {|\tilde{\mu}|} (\mathbb{R}^n) \le (1+\epsilon) \mathchoice {\left|\tilde{\mu}\right|} {|\tilde{\mu}|} {|\tilde{\mu}|} {|\tilde{\mu}|} (\mathbb{R}^n).\] Thus \(h\in A(\mathbb{R}^n)\) and \(\mathchoice {\left\lVert h\right\rVert} {\lVert h\rVert} {\lVert h\rVert} {\lVert h\rVert} _{W}\le (1+\epsilon) \mathchoice {\left\lVert\tilde{g}\right\rVert} {\lVert\tilde{g}\rVert} {\lVert\tilde{g}\rVert} {\lVert\tilde{g}\rVert} _{W}\). Letting \(\epsilon\to 0\) gives \[\label{eq:32wiener32norm32equiv} \mathchoice {\left\lVert g\right\rVert} {\lVert g\rVert} {\lVert g\rVert} {\lVert g\rVert} _{W(\Omega)} :=\;\inf_{\tilde{g}\in B(\mathbb{R}^n),\;\tilde{g}|_\Omega = g}\; \mathchoice {\left\lVert\tilde{g}\right\rVert} {\lVert\tilde{g}\rVert} {\lVert\tilde{g}\rVert} {\lVert\tilde{g}\rVert} _W =\;\inf_{\hat{h}\in C^{\infty}(\mathbb{R}^n),\;h\in A(\mathbb{R}^n),\;h|_\Omega = g}\; \mathchoice {\left\lVert h\right\rVert} {\lVert h\rVert} {\lVert h\rVert} {\lVert h\rVert} _W,\tag{74}\] completing the proof of 8 .

With not much effort, we can actually show the first infimum in 74 is attained:

9.3 Proof of 2↩︎

10 Details for the two-dimensional disc showcase↩︎

This appendix section contains the elementary calculations used in the disc showcase. The main text only uses the final coefficient formulas and the resulting second-order cone program.

10.1 Neumann Fourier–Bessel modes and the zero convention↩︎

We follow the DLMF convention for zeros of derivatives of Bessel functions [40]. Thus \(\beta_{p,\ell}\) is the \(\ell\)-th non-negative zero of \(J_p'\), \(\beta_{0,1}=0\), and all other \(\beta_{p,\ell}\) are positive.

The Bessel differential equation [40] implies \[\label{eq:32app32bessel32radial32eigen32equation} -\frac{1}{r}\frac{\,\mathrm{d}}{\,\mathrm{d}r}\left(r\frac{\,\mathrm{d}}{\,\mathrm{d}r} J_p(\beta r)\right)+\frac{p^2}{r^2}J_p(\beta r) =\beta^2J_p(\beta r).\tag{82}\] Therefore \[\Psi_{m,\ell}(r,\theta)=J_{|m|}(\beta_{|m|,\ell}r)e^{\mathrm{i}m\theta}\] is an eigenfunction of \(-\Delta\) on the unit disc. The boundary condition is Neumann because \[\partial_r J_{|m|}(\beta_{|m|,\ell}r)\big|_{r=1} =\beta_{|m|,\ell}J_{|m|}'(\beta_{|m|,\ell})=0.\]

The modes are orthogonal in \(L^2(D_1)\). The angular part gives \[\int_0^{2\pi}e^{\mathrm{i}m\theta}e^{-\mathrm{i}m'\theta}\,\mathrm{d}\theta =2\pi\delta_{m,m'}.\] The radial part is the usual Sturm–Liouville orthogonality for 82 . Thus, with \[N_{p,\ell}=\int_0^1J_p(\beta_{p,\ell}r)^2r\,\mathrm{d}r,\] we have \[\label{eq:32app32disc32FB32orthogonality} \int_0^{2\pi}\int_0^1 \Psi_{m,\ell}(r,\theta)\overline{\Psi_{m',\ell'}(r,\theta)}\,r\,\mathrm{d}r\,\mathrm{d}\theta =2\pi N_{|m|,\ell}\delta_{m,m'}\delta_{\ell,\ell'}.\tag{83}\] For the constant mode, \(N_{0,1}=1/2\). For \(\beta_{p,\ell}>0\), the closed form \[\label{eq:32app32disc32FB32norm32closed32form} N_{p,\ell} =\frac{1}{2}\left(1-\frac{p^2}{\beta_{p,\ell}^2}\right)J_p(\beta_{p,\ell})^2\tag{84}\] follows from the Bessel integral identity [40] after setting the upper limit equal to \(\beta_{p,\ell}\) and using \(J_p'(\beta_{p,\ell})=0\).

10.2 Jacobi–Anger expansion and the Wiener-type bound↩︎

The special-function input is the Jacobi–Anger expansion. It follows from the generating function [40]; the Jacobi–Anger expansions are listed explicitly in [40]. In the form used here, \[\label{eq:32app32JA32expansion32used} e^{-\mathrm{i}\rho r\cos(\theta-\varphi)} =\sum_{n\in\mathbb{Z}}(-\mathrm{i})^{|n|}J_{|n|}(\rho r)e^{\mathrm{i}n\theta}e^{-\mathrm{i}n\varphi}.\tag{85}\] Indeed, the more common version is \[e^{-\mathrm{i}z\cos\psi}=\sum_{n\in\mathbb{Z}}(-\mathrm{i})^nJ_n(z)e^{\mathrm{i}n\psi},\] and 85 follows by taking \(\psi=\theta-\varphi\) and using \(J_{-p}=(-1)^pJ_p\) for integer \(p\) [40].

Multiplying 85 by \(e^{\mathrm{i}m\varphi}\) and integrating in \(\varphi\) gives \[\begin{align} \int_0^{2\pi}e^{-\mathrm{i}\rho r\cos(\theta-\varphi)}e^{\mathrm{i}m\varphi}\,\mathrm{d}\varphi &=\sum_{n\in\mathbb{Z}}(-\mathrm{i})^{|n|}J_{|n|}(\rho r)e^{\mathrm{i}n\theta} \int_0^{2\pi}e^{-\mathrm{i}n\varphi}e^{\mathrm{i}m\varphi}\,\mathrm{d}\varphi \\ &=2\pi(-\mathrm{i})^{|m|}J_{|m|}(\rho r)e^{\mathrm{i}m\theta}. \end{align}\] Hence \[\label{eq:32app32disc32FB32plane32wave32average} J_{|m|}(\rho r)e^{\mathrm{i}m\theta} =\frac{\mathrm{i}^{|m|}}{2\pi} \int_0^{2\pi} e^{-\mathrm{i}\rho(x\cos\varphi+y\sin\varphi)}e^{\mathrm{i}m\varphi}\,\mathrm{d}\varphi.\tag{86}\] Thus each Fourier–Bessel mode is an average of plane waves with frequencies on the circle \(|\boldsymbol{\xi}|=\rho\). The factor \(\mathrm{i}^{|m|}\) has modulus one, and the averaging measure in 86 has total variation one. Therefore, in the Fourier–Stieltjes sense explained in 9.1, \[\mathchoice {\left\lVert\Psi_{m,\ell}\right\rVert} {\lVert\Psi_{m,\ell}\rVert} {\lVert\Psi_{m,\ell}\rVert} {\lVert\Psi_{m,\ell}\rVert} _{W(D_1)}\le 1.\]

10.3 The product integral and the plane-wave coefficient matrix↩︎

In this subsection, we give a closed-form formula for the coefficient of each plane wave in the Fourier–Bessel basis.

For \((j,k)\in\mathcal{G}_R\), write \[\omega_{j,k}=\frac{\pi}{2}(j,k), \qquad q_{j,k}=|\omega_{j,k}|, \qquad \varphi_{j,k}=\mathop{\mathrm{Arg}}(j+\mathrm{i}k).\] For \((j,k)=(0,0)\), the value of \(\varphi_{0,0}\) is arbitrary. Then \[e^{-\mathrm{i}\omega_{j,k}\cdot(x,y)} =e^{-\mathrm{i}q_{j,k}r\cos(\theta-\varphi_{j,k})}.\] We calculate the Fourier–Bessel coefficients by the inner products \[\label{eq:32disc32FB32residual32coeff32P} P_{m,\ell}^{j,k} =\frac{1}{2\pi N_{|m|,\ell}} \int_0^{2\pi}\int_0^1 e^{-\mathrm{i}q_{j,k}r\cos(\theta-\varphi_{j,k})} J_{|m|}(\beta_{|m|,\ell}r)e^{-\mathrm{i}m\theta}\,r\,\mathrm{d}r\,\mathrm{d}\theta.\tag{87}\] Substituting the Jacobi–Anger expansion 85 into the coefficient definition 87 , the angular integral becomes \[\frac{1}{2\pi}\int_0^{2\pi}e^{\mathrm{i}n\theta}e^{-\mathrm{i}m\theta}\,\mathrm{d}\theta =\delta_{n,m},\] so only the term \(n=m\) remains. The remaining radial integral gives \[P_{m,\ell}^{j,k} =(-\mathrm{i})^{|m|}e^{-\mathrm{i}m\varphi_{j,k}} \frac{1}{N_{|m|,\ell}}\int_0^1 J_{|m|}(q_{j,k}r) J_{|m|}(\beta_{|m|,\ell}r) r \,\mathrm{d}r.\] The radial integral can be simplified using the standard Bessel product identity [40]: \[\label{eq:32app32disc32FB32H32closed32form} \int_0^1J_p(qr)J_p(\beta r)r\,\mathrm{d}r =\frac{qJ_{p+1}(q)J_p(\beta)-\beta J_p(q)J_{p+1}(\beta)}{q^2-\beta^2}.\tag{88}\] At \(q=\beta\), the right-hand side is understood by continuity.

10.4 A practical extension bound for the local Wiener norm on disc↩︎

The Fourier–Bessel coefficient certificate above is straightforward to incorporate into the SOCP. However, it is computationally expensive because of its slow polynomial convergence. For a fixed residual, a direct Sobolev-type extension bound can give a much less expensive a posteriori estimate of the local Wiener norm, and can sometimes be sharper than the Fourier–Bessel coefficient certificate. We record the version used in the numerical checks. It is independent of the Fourier–Bessel basis.

Let \(g\) be a smooth function on \(D_1\), and let \(a>0\) be a free scaling parameter. By the same Cauchy–Schwarz argument as in 9 , but with the weight \(1+a^2 \mathchoice {\left|\boldsymbol{\xi}\right|} {|\boldsymbol{\xi}|} {|\boldsymbol{\xi}|} {|\boldsymbol{\xi}|} ^2\), every extension \(G\in H^2(\mathbb{R}^2)\) of \(g\) satisfies \[\label{eq:32app32practical32scaled32Sobolev32Wiener} \mathchoice {\left\lVert g\right\rVert} {\lVert g\rVert} {\lVert g\rVert} {\lVert g\rVert} _{W(D_1)} \le \mathchoice {\left\lVert G\right\rVert} {\lVert G\rVert} {\lVert G\rVert} {\lVert G\rVert} _{W(\mathbb{R}^2)} \le \frac{1}{2a\sqrt{\pi}} \mathchoice {\left\lVert G-a^2\Delta G\right\rVert} {\lVert G-a^2\Delta G\rVert} {\lVert G-a^2\Delta G\rVert} {\lVert G-a^2\Delta G\rVert} _{L^2(\mathbb{R}^2)} .\tag{89}\] Indeed, \[\int_{\mathbb{R}^2}\frac{\,\mathrm{d}\boldsymbol{\xi}}{(1+a^2 \mathchoice {\left|\boldsymbol{\xi}\right|} {|\boldsymbol{\xi}|} {|\boldsymbol{\xi}|} {|\boldsymbol{\xi}|} ^2)^2}=\frac{\pi}{a^2},\] and Plancherel’s theorem gives \[\mathchoice {\left\lVert(1+a^2 \mathchoice {\left|\boldsymbol{\xi}\right|} {|\boldsymbol{\xi}|} {|\boldsymbol{\xi}|} {|\boldsymbol{\xi}|} ^2)\widehat G\right\rVert} {\lVert(1+a^2 \mathchoice {\left|\boldsymbol{\xi}\right|} {|\boldsymbol{\xi}|} {|\boldsymbol{\xi}|} {|\boldsymbol{\xi}|} ^2)\widehat G\rVert} {\lVert(1+a^2 \mathchoice {\left|\boldsymbol{\xi}\right|} {|\boldsymbol{\xi}|} {|\boldsymbol{\xi}|} {|\boldsymbol{\xi}|} ^2)\widehat G\rVert} {\lVert(1+a^2 \mathchoice {\left|\boldsymbol{\xi}\right|} {|\boldsymbol{\xi}|} {|\boldsymbol{\xi}|} {|\boldsymbol{\xi}|} ^2)\widehat G\rVert} _{L^2(\mathbb{R}^2)} = \mathchoice {\left\lVert G-a^2\Delta G\right\rVert} {\lVert G-a^2\Delta G\rVert} {\lVert G-a^2\Delta G\rVert} {\lVert G-a^2\Delta G\rVert} _{L^2(\mathbb{R}^2)} .\] Thus the problem reduces to choosing an extension that makes the last norm in 89 small.

We now compute the optimal exterior contribution on the unit disc. Write the boundary trace and the radial derivative of \(g\) as \[\label{eq:32app32practical32boundary32Fourier32data} g(1,\theta)=\sum_{n\in\mathbb{Z}}c_ne^{\mathrm{i}n\theta}, \qquad \partial_rg(1,\theta)=\sum_{n\in\mathbb{Z}}d_ne^{\mathrm{i}n\theta}.\tag{90}\] The derivative in 90 is the outward radial derivative from the disc. Any \(H^2\) extension that equals \(g\) on \(D_1\) must have the same trace and radial derivative at \(r=1\).

Let \(E=\mathbb{R}^2\setminus \overline{D_1}\) and minimize \[J_a(u)=\int_E \mathchoice {\left|u-a^2\Delta u\right|} {|u-a^2\Delta u|} {|u-a^2\Delta u|} {|u-a^2\Delta u|} ^2\,\mathrm{d}x\] among decaying exterior functions with the boundary data 90 . The Euler–Lagrange equation is \[\label{eq:32app32practical32EL32equation} (1-a^2\Delta)^2u=0\qquad \text{in }E.\tag{91}\] Set \(\kappa=1/a\). The decaying separated solutions of the modified Helmholtz equation are \(K_{ \mathchoice {\left|n\right|} {|n|} {|n|} {|n|} }(\kappa r)e^{\mathrm{i}n\theta}\), where \(K_p\) is the modified Bessel function of the second kind; see the modified Bessel equation and its standard solutions in [40]. Since 91 contains the square of the modified Helmholtz operator, a second decaying radial solution is obtained by differentiating \(K_{ \mathchoice {\left|n\right|} {|n|} {|n|} {|n|} }(\kappa r)\) with respect to \(\kappa\), namely \(rK_{ \mathchoice {\left|n\right|} {|n|} {|n|} {|n|} }'(\kappa r)\). Hence \[\label{eq:32app32practical32exterior32solution} u(r,\theta) =\sum_{n\in\mathbb{Z}}\left(A_nK_{ \mathchoice {\left|n\right|} {|n|} {|n|} {|n|} }(\kappa r) +B_n rK_{ \mathchoice {\left|n\right|} {|n|} {|n|} {|n|} }'(\kappa r)\right)e^{\mathrm{i}n\theta}, \qquad r>1.\tag{92}\] Here and below, the prime denotes differentiation with respect to the argument of \(K_{ \mathchoice {\left|n\right|} {|n|} {|n|} {|n|} }\).

For brevity set \[K=K_{ \mathchoice {\left|n\right|} {|n|} {|n|} {|n|} }(\kappa), \qquad K'=K_{ \mathchoice {\left|n\right|} {|n|} {|n|} {|n|} }'(\kappa), \qquad u_n(\kappa)=-\frac{K'}{K}.\] Using the modified Bessel equation at \(\kappa\), \[\kappa^2K''+\kappa K'-(\kappa^2+n^2)K=0,\] the boundary conditions \(u(1,\theta)=g(1,\theta)\) and \(\partial_ru(1,\theta)=\partial_rg(1,\theta)\) give, mode by mode, \[\label{eq:32app32practical32AB32system} A_nK+B_nK'=c_n, \qquad \kappa A_nK'+B_n\frac{\kappa^2+n^2}{\kappa}K=d_n.\tag{93}\] Solving the two-by-two system gives \[\label{eq:32app32practical32Bn32formula} B_n=-\frac{\kappa\left(d_n+\kappa u_n(\kappa)c_n\right)}{K\left(\kappa^2u_n(\kappa)^2-\kappa^2-n^2\right)} .\tag{94}\] The nontrivial factor in the denominator is positive: \[\label{eq:32app32practical32positive32denominator} \kappa^2u_n(\kappa)^2-\kappa^2-n^2>0.\tag{95}\] This follows from the integral identity [40] \[\label{eq:32app32practical32K32square32integral} \int_\kappa^\infty K_{ \mathchoice {\left|n\right|} {|n|} {|n|} {|n|} }(t)^2t\,\mathrm{d}t =\frac{\kappa^2}{2} \left(K_{ \mathchoice {\left|n\right|} {|n|} {|n|} {|n|} }'(\kappa)^2- \left(1+\frac{n^2}{\kappa^2}\right)K_{ \mathchoice {\left|n\right|} {|n|} {|n|} {|n|} }(\kappa)^2\right).\tag{96}\]

Let \[v=u-a^2\Delta u.\] The first term in 92 is annihilated by \(1-a^2\Delta\). Differentiating the modified Helmholtz equation with respect to \(\kappa\) shows that \[(1-a^2\Delta)\big(rK_{ \mathchoice {\left|n\right|} {|n|} {|n|} {|n|} }'(\kappa r)e^{\mathrm{i}n\theta}\big) =-2aK_{ \mathchoice {\left|n\right|} {|n|} {|n|} {|n|} }(\kappa r)e^{\mathrm{i}n\theta}.\] Therefore \[v(r,\theta) =-2a\sum_{n\in\mathbb{Z}}B_nK_{ \mathchoice {\left|n\right|} {|n|} {|n|} {|n|} }(\kappa r)e^{\mathrm{i}n\theta}.\] Using orthogonality in \(\theta\), the change of variables \(t=\kappa r\), the integral identity 96 , and the expression 94 , we obtain \[\begin{align} \label{eq:32app32practical32exterior32cost} \inf_u J_a(u) &= \mathchoice {\left\lVert v\right\rVert} {\lVert v\rVert} {\lVert v\rVert} {\lVert v\rVert} _{L^2(E)}^2 \notag\\ &=2\pi a^2\sum_{n\in\mathbb{Z}} \frac{2 \mathchoice {\left|u_n(1/a)c_n+a d_n\right|} {|u_n(1/a)c_n+a d_n|} {|u_n(1/a)c_n+a d_n|} {|u_n(1/a)c_n+a d_n|} ^2}{u_n(1/a)^2-(na)^2-1}. \end{align}\tag{97}\] Combining the interior contribution on \(D_1\) with the optimal exterior contribution, 89 gives the practical bound \[\label{eq:32app32practical32Wiener32final} \mathchoice {\left\lVert g\right\rVert} {\lVert g\rVert} {\lVert g\rVert} {\lVert g\rVert} _{W(D_1)} \le \frac{1}{2a\sqrt{\pi}} \left( \mathchoice {\left\lVert g-a^2\Delta g\right\rVert} {\lVert g-a^2\Delta g\rVert} {\lVert g-a^2\Delta g\rVert} {\lVert g-a^2\Delta g\rVert} _{L^2(D_1)}^2 +2\pi a^2\sum_{n\in\mathbb{Z}} \frac{2 \mathchoice {\left|u_n(1/a)c_n+a d_n\right|} {|u_n(1/a)c_n+a d_n|} {|u_n(1/a)c_n+a d_n|} {|u_n(1/a)c_n+a d_n|} ^2}{u_n(1/a)^2-(na)^2-1} \right)^{1/2} .\tag{98}\] The parameter \(a>0\) is free and can be optimized numerically. The case \(a=1\) is the unscaled bound obtained from the operator \(1-\Delta\).

References↩︎

---

1. hongkang@stanford.edu↩︎

2. lexing@stanford.edu↩︎

3. This work is partially supported by the U.S. Department of Energy, Office of Science, Accelerated Research in Quantum Computing Centers, Quantum Utility through Advanced Computational Quantum Algorithms, grant no. DE-SC0025572.↩︎