Getting started

Installation

felcu needs Python 3.13+ and NumPy 2.2+. Install from a clone of the repository:

pip install .

For development (tests, plus an editable install so changes take effect immediately):

pip install -e ".[dev]"

A first fit

Fit \(e^{\mathcal H}\) for a Hermitian \(\mathcal H\) with \(\operatorname{spec}(\mathcal H) \subseteq [0.1, 2]\) to a target tolerance \(\epsilon_{\rm t} = 10^{-10}\):

from felcu import solve_fe_coefficients

result = solve_fe_coefficients(
    "exponentiation",
    sigma_min=0.1,
    sigma_max=2.0,
    target_error=1e-10,
)

print(result.converged)            # True if target_error was met
print(result.m)                    # truncation order chosen
print(result.achieved_linf_error)  # epsilon_m on the fitted set
print(result.alpha)                # subnormalisation alpha_m

The solver tries a range of truncation orders \(m\) and, among those that meet the target, returns the one with the smallest subnormalisation \(\alpha_m\). Pass selection="smallest_m" for the shortest series instead. Every candidate it tried is kept in result.all_results.

If no candidate meets the target, the solver still returns the best one, with converged=False, and issues a felcu.ConvergenceWarning. Widen m_values or loosen target_error, or silence the warning with warnings.simplefilter("ignore", felcu.ConvergenceWarning) when missing the target is expected.

Choosing the range of m

By default the solver sweeps \(m\) up to 256. Because \(\alpha\) usually keeps falling as \(m\) grows, the default smallest_alpha selection often picks the top of the range. If circuit depth matters, cap the range, e.g. m_values=range(4, 41, 2). See Tuning the solver.

This uses the default Sobolev-regularised fit. Pass method="reflected" or method="arcsine_taylor" for the paper's explicit extensions, which trade convergence against subnormalisation differently; see Method.

Applying the result to a matrix

The fit is done in a scaled variable \(x\), so the operator you implement is the series \(g_m\) evaluated at the scaled matrix \(\mathcal G\):

\[ f(\mathcal H) \approx g_m(\mathcal G), \qquad \mathcal G = \tau(\mathcal H - \mu I), \]

with \(\mu\) the centre of the spectrum (result.mu) and \(\tau\) the scale (result.tau). For singular values, \(\mathcal H\) is the Hermitian dilation of \(A\) and \(\mu = 0\). The result applies the scaling for you, so you can check the approximation classically:

from felcu import solve_fe_coefficients

H = ...  # Hermitian matrix
result = solve_fe_coefficients("exponentiation", matrix=H, target_error=1e-8)

result.evaluate(x)   # ~= exp(x) for x in the spectral interval
result.apply(H)      # ~= expm(H), by eigendecomposition

apply needs a Hermitian matrix; a non-Hermitian one is block encoded through its Hermitian dilation, so apply the result to that. To evaluate \(g_m\) itself in the scaled variable, use evaluate_fourier_series.

Passing matrix= instead of sigma_min/sigma_max lets felcu read the spectral interval from the matrix itself. One or the other is required.

Coefficients for an LCU circuit

The real coefficients a0, ak, bk describe \(g_m\) in the sine/cosine basis. For a circuit you usually want the exponential form

\[ g_m(x) = \sum_{k=-m}^{m} c_k e^{ikx}, \qquad c_0 = a_0,\quad c_{\pm k} = \frac{b_k \mp i a_k}{2}, \]

which is available as result.complex_coefficients, ordered \([c_{-m}, \dots, c_0, \dots, c_m]\). You can also convert any (a0, ak, bk) triple with to_complex_coefficients.

This order, by \(k\), is the one to use for evaluating the series, and the one saved by to_json. The compressed circuit takes the same coefficients in a different order, with a zero after \(c_0\):

\[ [c_0,\ 0,\ c_1,\ c_{-1},\ \dots,\ c_m,\ c_{-m}], \]

\(2m + 2\) entries in all. The zero is required by the compressed circuit and must be kept. This order is result.circuit_coefficients, and it is what to pass to the circuit, not complex_coefficients. Circuit code that imports felcu can take the array directly:

from felcu import solve_fe_coefficients

coefficients = solve_fe_coefficients("exponentiation", sigma_min=0.1, sigma_max=2.0).circuit_coefficients

To paste the coefficients somewhere by hand, result.format_circuit_coefficients() gives them as a bare Python list, zero included:

print(result.format_circuit_coefficients())

[(3.3903633881883226+0j), 0j, (-0.06034330806667995-1.1020475420922582j), (-0.06034330806667995+1.1020475420922582j), ...]

Saving and loading results

A result can be written to JSON and read back exactly:

from felcu import FourierExtensionResult

result.to_json("exp_coefficients.json")
loaded = FourierExtensionResult.from_json("exp_coefficients.json")

The file is standard JSON, so it can also be read outside Python. It holds every field of the result: the parameters (problem, method, eta, tau, sobolev_r, ...), the errors and \(\alpha_m\), the real coefficients a0, ak, bk, and the complex coefficients as

"complex_coefficients": {"k": [-m, ..., m], "real": [...], "imag": [...]}

Floats are written with enough digits to round-trip exactly; an infinite \(\tau\) (a zero-width interval) is stored as the string "inf", and the Sobolev-only fields (lambda_used, sobolev_r, sobolev_w) are null for the explicit methods. Pass include_candidates=True to also keep all_results, every order tried in the sweep. to_dict and from_dict give the same data as a Python dictionary.