Fourier LCU for Non-Unitary Decompositions

• The paper introduces Fourier-LCU, a framework that decomposes non-unitary operators into linear combinations of unitaries using periodic extension and Fourier sine series for exponential error convergence.
• It converts sine series into complex-exponential forms to map operator terms onto pairs of unitaries, facilitating a block-encoding method with double-logarithmic subnormalization scaling.
• The approach employs convex optimization for coefficient regularization, achieving a Pareto-optimal trade-off between error tolerance and resource efficiency in quantum algorithm implementations.

A Fourier Linear Combination of Unitaries (Fourier-LCU) is a general analytic method for decomposing arbitrary non-unitary operators into accurate, exponentially convergent linear combinations of unitary operators. This is accomplished via smooth periodic extension and Fourier sine series techniques, yielding a block-encoding whose subnormalization parameter exhibits double-logarithmic scaling in the target error. The framework leverages convex optimization to regularize the coefficients for specific error budgets, tracing out a Pareto front for subnormalization-versus-error. These advances constitute a versatile approach for non-unitary quantum algorithms and circuits (Brearley et al., 25 Jan 2026).

1. Periodic Extension and Fourier Sine Series Construction

To represent an operator via a LCU, the core technical step is constructing a periodic extension of the identity function $f(\tau) = \tau$ within a given interval. Fixing $\tau \in [-\pi/\eta, \pi/\eta]$ for some $\eta > 1$, one extends $f(\tau)$ to $[−\pi, \pi]$ as a $2\pi$-periodic, odd, and infinitely differentiable function. This ensures analyticity and supports exponential Fourier coefficient decay.

The extension yields a truncated $m$-term sine series approximation:

$f(\tau) \approx \sum_{k=1}^{m} a_k \sin(k \tau)$

The optimal coefficients $a_k$ are determined by a continuous least-squares problem over $[−\pi/\eta, \pi/\eta]$:

$$a^{\rm LS} = \arg\min_{a \in \mathbb{R}^m}\int_{-\pi/\eta}^{\pi/\eta} \left( \tau - \sum_{k=1}^m \tilde{a}_k \sin(k\tau) \right)^2 d\tau$$

The resulting normal equations are $G(\eta) a = b(\eta)$, with

$$G_{jk} = \int_{-\pi/\eta}^{\pi/\eta} \sin(j\tau)\sin(k\tau)\,d\tau, \quad b_j = \int_{-\pi/\eta}^{\pi/\eta} \tau\,\sin(j\tau)\,d\tau$$

Since $\eta > 1$, the coefficients $a_k$ exhibit exponential decay in $k$, ensuring rapid convergence.

2. Complex-Exponential Formulation and Unitary Mapping

The sine series can be rewritten using the Euler identity:

$\sin(kx) = \frac{e^{ikx} - e^{-ikx}}{2i}$

Substitution yields:

$$f(x) = \sum_{k=1}^m a_k \sin(kx) = \sum_{n=-m}^{m} c_n e^{inx}$$

where $c_{+k} = -\frac{i}{2} a_k$, $c_{-k} = +\frac{i}{2} a_k$, and $c_0 = 0$. Consequently, each sine term maps to a pair of unitaries with complex weights, forming the desired LCU structure.

3. Application to Arbitrary Non-Unitary Operators

Let $A$ be a general (potentially non-unitary) operator. Decompose $A$ into Hermitian and anti-Hermitian components:

$$A = H_1 + i H_2,\quad H_1 = \frac{A + A^\dagger}{2},\quad H_2 = \frac{A - A^\dagger}{2i}$$

The LCU approximation proceeds by:

• Choosing $\tau$ so that $$\text{spec}(\tau H_1) \subseteq [−\pi/\eta, \pi/\eta]$$;
• Approximating $H_1$ and $iH_2$ by sine series as above;
• Rewriting each sine term in unitary difference form.

The full LCU for $A$ to exponential error $O(e^{-c m})$ is:

$$A \approx \sum_{k=1}^m \frac{a_k}{2\tau} \left[ i e^{-i k \tau H_1} - i e^{+i k \tau H_1} + e^{+i k \tau H_2} - e^{-i k \tau H_2} \right] = \sum_{j=1}^{4m} \kappa_j U_j$$

Each $U_j$ is one of $e^{\pm i k \tau H_1}, e^{\pm i k \tau H_2}$ with real weights $\kappa_j$ given by $\pm a_k/(2\tau)$ (up to phase).

4. Block-Encoding and Subnormalization Scaling

Employing standard LCU block-encoding [Childs–Wiebe 2012], one introduces $n_a = \lceil \log_2(4m) \rceil$ ancillas, prepares amplitude state $V$, applies controlled-$U_j$ gates, and uncomputes via $W^\dagger$:

$$(\langle 0|^{\otimes n_a} \otimes I) [ W^\dagger (\mathrm{ctrl}\text{-}U) V ] (|0\rangle^{\otimes n_a} \otimes I) = A/\alpha + O(\epsilon)$$

with subnormalization parameter

$$\alpha = \sum_{j=1}^{4m} | \kappa_j | = \frac{2}{\tau} \sum_{k=1}^{m} | a_k | = \frac{2\eta}{\pi} \max(\|H_1\|_2,\|H_2\|_2) \sum_{k=1}^{m} | a_k |$$

Since $a_k$ decay exponentially and empirically $\sum_{k=1}^{m} |a_k| = O(\log m)$, the total normalization satisfies:

$\alpha = O(\|A\|_2 \cdot \log \log(1/\epsilon))$

This double-logarithmic scaling in $1/\epsilon$ is a substantial improvement over previous polynomial relationships between subnormalization and error.

5. Coefficient Regularization and Pareto Front Optimization

Because the sine dictionary is overcomplete for $\eta > 1$, there exist infinitely many coefficient sets yielding nearly identical $L^2$ error yet different $L^1$ summations (impacting $\alpha$). Regularization is performed via convex optimization, exploiting the trade-off:

$$J(a; \lambda) = \|\tau - \Phi a\|_2 + \lambda \frac{2\eta}{\pi} \|a\|_1$$

for $a \in \mathbb{R}^m$ and dictionary matrix $\Phi_{jk} = \sin(k\tau_j)$. At fixed error budget $\epsilon$, the $L^1$-minimization is:

$$\alpha^*(\epsilon) = \min\left\{ \frac{2\eta}{\pi} \|a\|_1 : \|\tau - \Phi a\|_2 \leq \epsilon \right\}$$

Standard convex solvers or homotopy/LASSO-type path tracking yield the unique Pareto front $\alpha(\epsilon)$. It is proven that $J^*_m(\lambda)$ and $\alpha^*_m(\epsilon)$ are nonincreasing in $m$ and converge to a finite limit as $m \to \infty$. Numerically, “sweeping” $\lambda$ down to zero identifies the lowest possible $\alpha$ at target $\epsilon$.

6. Implementation Procedures and Practical Implications

The Fourier-LCU methodology is summarized by the following stepwise procedure:

1. Construct an analytic $2\pi$-periodic, odd extension of $f(\tau) = \tau$ on $[−\pi/\eta, \pi/\eta]$.
2. Compute the exponentially convergent truncated sine series via least squares.
3. Convert each $\sin(k \tau H_i)$ into LCUs of $e^{\pm i k \tau H_i}$ for $i = 1,2$; assemble $A = H_1 + i H_2$ via weighted sums.
4. Realize the decomposition as an $(\alpha, n_a, \epsilon)$ block-encoding, utilizing $\alpha$ scaling as $O(\log \log(1/\epsilon))$.
5. Optionally, re-optimize coefficients with $\ell_1$-regularized least squares to minimize $\alpha$ at fixed $\epsilon$, thereby mapping out the Pareto front. All essential equations for $a_k$, $c_n$, $\kappa_j$, and $\alpha$ as well as the regularization trade-offs are explicitly stated; optimized $a_k$ values for various $m$ are tabulated in Table B of the corresponding source (Brearley et al., 25 Jan 2026).

A plausible implication is that non-unitary quantum algorithms leveraging Fourier-LCU can reach error targets at far lower resource cost than via polynomial-scaling frameworks, and coefficient regularization can yield trainable sparsity for practical block-encodings.