Coefficients#

A characteristic property of any Fourier model are its coefficients. Our package can, given a model, calculate the corresponding coefficients.

In the simplest case, this could look as follows:

from qml_essentials.model import Model
from qml_essentials.coefficients import Coefficients

model = Model(
            n_qubits=2,
            n_layers=1,
            circuit_type="Hardware_Efficient",
        )

coeffs, freqs = Coefficients.get_spectrum(model)

Here, the coefficients are stored in the coeffs variable, and the corresponding frequency indices are stored in the freqs variable.

But wait! There is much more to this. Let's keep on reading if you're curious .

Detailled Explanation#

To visualize what happens, let's create a very simplified Fourier model

class Model_Fct:
    def __init__(self, c, f):
        self.c = c
        self.f = f
        self.degree = (2*max(f)+1,)
        self.frequencies = f
        self.n_input_feat = 1

    def __call__(self, inputs, **kwargs):
        return np.sum([c * np.exp(-1j * inputs * f) for f, c in zip(self.f, self.c)], axis=0)

This model takes a vector of coefficients and frequencies on instantiation. When called, these coefficients and frequencies are used to compute the output of the model, which is the sum of sine functions determined by the length of the vectors. Let's try that for just two frequencies:

freqs = [-3, -1.5, 0, 1.5, 3]
coeffs = [1, 1, 0, 1, 1]

fs = max(freqs) * 2 + 1
model_fct = Model_Fct(coeffs, freqs)

x = np.arange(0, 2 * np.pi, 2 * np.pi / fs)
out = model_fct(x)

We can now calculate the Fast Fourier Transform of our model:

X = np.fft.fft(out) / len(out)
X_shift = np.fft.fftshift(X)
X_freq = np.fft.fftfreq(X.size, 1/fs)
X_freq_shift = np.fft.fftshift(X_freq)

You may notice, that something isn't quite right here; we specified the frequencies [1.5,3] earlier, but get frequencies for [0,1,2,3]. This is because, we chose the wrong resolution for the FFT, i.e. the window length was too short. In our framework we can achieve the same and above, while simultanously applying the fix to our problem, i.e. setting mts=2. This additional variable effectively doubles the window length which gives us then the possibility to obtain frequencies "between" the integer valued frequencies seen above.

X_shift, X_freq_shift = Coefficients.get_spectrum(model_fct, mts=2, shift=True)

Note, that applying the shift can be controlled with the optional shift argument.

Another important point is, that the force_mean flag is set, and the execution_type is is implicitly set to expval. This is mainly because, we require a single expectation value to calculate the coefficients.

Increasing the Resolution#

You might have noticed that we choose our sampling frequency fs in such a way, that it just fulfills the Nyquist criterium. Also the number of samples x are just enough to sufficiently represent our function. In such a simplified scenario, this is fine, but there are cases, where we want to have more information both in the time and frequency domain. Therefore, two additional arguments exist in the get_spectrum method: - mfs: The multiplier for the highest frequency. Increasing this will increase the width of the spectrum - mts: The multiplier for the number of time samples. Increasing this will increase the resolution of the time domain and therefore "add" frequencies in between our original frequencies. - trim: Whether to remove the Nyquist frequency if spectrum is even. This will result in a symmetric spectrum

X_shift, X_freq_shift = Coefficients.get_spectrum(model_fct, mfs=2, mts=3, shift=True)

Note that, as the frequencies change with the mts argument, we have to take that into account when calculating the frequencies with the last call.

Feel free to checkout our jupyter notebook if you would like to play around with this.

A sidenote on the performance; Increasing the mts value effectively increases the input lenght that goes into the model. This means that mts=2 will require twice the time to compute, which will be very noticable when running noisy simulations.

Power spectral density#

In some cases it can be useful to get the power spectral density (PSD). As calculation of this metric might differ between the different research domains, we included a function to get the PSD of a given spectrum using the following formula:

where \(F\) is the spectrum and \(n_\text{samples}\) the length of the input vector.

model = Model(
    n_qubits=4,
    n_layers=1,
    circuit_type="Circuit_19",
    random_seed=1000
)

coeffs, freqs = Coefficients.get_spectrum(model, mfs=1, mts=1, shift=True)

psd = Coefficients.get_psd(coeffs)

Evaluating a Fourier Series#

Given a set of coefficients and frequencies (e.g. those returned by get_spectrum), the corresponding Fourier series can be evaluated at arbitrary input points with evaluate_Fourier_series:

coeffs, freqs = Coefficients.get_spectrum(model)

y = Coefficients.evaluate_Fourier_series(coeffs, freqs, inputs=0.5)

This accepts both single- and multi-dimensional frequency arrays, mirroring the shapes produced by get_spectrum.

Analytic Coefficients#

All of the calculations above were performed by applying a Fast Fourier Transform to the output of our Model. However, we can also calculate the coefficients analytically.

This can be achieved by the so called FourierTree class:

from qml_essentials.coefficients import FourierTree

fourier_tree = FourierTree(model)
an_coeffs, an_freqs = fourier_tree.get_spectrum(force_mean=True)

Unlike the FFT, this gives us the precise coefficients, solely depending on the parameters. We can verify this by comparing it to the previous results:

Technical Details#

We use an approach developed by Nemkov et al., which was later extended by Wiedmann et al.. The implementation is also inspired by the corresponding code for Nemkov et al.'s paper.

In Nemkov et al.'s algorithm the first step is to separate Clifford and non-Clifford gates, such that all Clifford gates can be regarded as part of the observable, and the actual circuit only consists of Pauli rotations (cf. qml_essentials.pauli.PauliCircuit). The main idea is then to split each Pauli rotation into sine and cosine product terms to obtain the coefficients, which are only dependent on the parameters of the circuit.

The Clifford commutation and Pauli bookkeeping are implemented symbolically on a stabilizer-tableau Pauli representation (jaqsi.paulis.PauliWord), and the parameter-dependent coefficients are evaluated with vectorized operations. This makes the tree fast to build and evaluate, and it supports multiple input features: in that case get_spectrum returns, per observable, the multi-dimensional frequency vectors (shape (n_freqs, n_features)) and their coefficients.

The symbolic core is gate-set agnostic and lives in jaqsi.paulis:

• PauliWord — an n-qubit Pauli in the symplectic \(i^{\text{phase}}\,X^x Z^z\) representation, with compose, commutes_with, conjugate_by_clifford (Clifford tableau evolution), and a to_matrix/from_matrix bridge to dense operators.
• Operation.is_clifford — a class flag marking the standard Clifford gates (\(I, X, Y, Z, H, S, \text{CX}, \text{CY}, \text{CZ}, \text{SWAP}\)); conjugate_by_clifford uses fast symbolic rules for the common ones and an exact matrix fallback for the rest.
• Operation.decompose() — expresses composite gates (Rot, CRX/CRY/CRZ, CZ) in terms of Clifford + Pauli-rotation primitives. Diagonal-Hamiltonian encodings (DiagonalQubitUnitary, e.g. the Golomb strategy) are expanded into commuting Pauli-\(Z\) rotations via the Walsh-Hadamard transform of their real generator.

Custom Circuits and Input Scaling#

The FourierTree operates on whatever circuit the model's script records, so a custom variational circuit can be analysed by replacing model.script:

import jaqsi as js
from jaqsi import Gates as g

def variational(params, inputs, *args, **kwargs):
    params = params.squeeze()
    g.PauliRot(1 * inputs, "Y", wires=[0])
    g.PauliRot(params[0], "XZX", wires=[0, 1, 2])
    g.PauliRot(3 * inputs, "XZ", wires=[0, 1])
    g.PauliRot(params[1], "YY", wires=[0, 1])
    g.PauliRot(9 * inputs, "XY", wires=[0, 1])

model = Model(n_qubits=3, n_layers=1, observables=0)
model._params_shape = (2, 1)
model.initialize_params()       # reallocate params for the custom shape
model.script = js.Script(f=variational, n_qubits=3)

coeffs, freqs = FourierTree(model).get_spectrum()

The tree distinguishes encoding from variational rotations automatically. Because every canonical rotation angle is an affine function of the inputs (encodings apply \(\omega_k\,x_f\); Clifford commutation only flips a generator's sign), the tree perturbs each input feature in turn and reads off, from the change in the recorded angles, which rotations depend on which feature and with what signed integer scaling \(\omega_k\).

The only requirement is that each encoding rotation be linear in a single feature, \(\omega_k\,x_f\). Heterogeneous scalings such as \((1, 3, 9)\) are then resolved to the correctly scaled frequencies (here \(\{-13, -11, -7, -5, 5, 7, 11, 13\}\)) instead of unit counts. Scalings may be rational rather than integer for some encodings: a diagonal-Hamiltonian encoding \(\exp(-i x\,\mathrm{diag}(h))\) (such as the Golomb strategy) decomposes into commuting Pauli-\(Z\) rotations whose scalings are dyadic rationals, while the resulting model frequencies stay integer (the pairwise differences of the Hamiltonian eigenvalues). A rotation depending on more than one feature raises NotImplementedError; and the method="dp" exact-support path does not support non-unit scaling (it raises for Golomb encodings — use method="tree").

For the numerical FFT (Coefficients.get_spectrum) a custom circuit additionally requires model.degree to be set high enough to resolve the largest frequency. The sampling grid is built from model.degree (the number of frequencies, i.e. \(2\,\omega_\text{max} + 1\)), not model.frequencies. Because model.degree reflects the nominal integer comb and not any scaling applied through enc_params or enc_pulse_params, a model whose effective frequencies exceed the nominal degree by a factor \(s > 1\) needs mfs at least \(\lceil s \rceil\) to avoid aliasing the scaled components.

Estimating the Exact Spectrum#

The number of frequencies a model can represent is, by default, estimated naively from the encoding (see model.frequencies / model.degree). This estimate is an upper bound: as shown by Wiedmann et al., some coefficients are constrained to zero for all parameter values, so the true spectrum can be smaller.

The FourierTree lets us obtain the exact spectrum of a model via the opt-in Model.exact_spectrum method:

model = Model(
    n_qubits=3,
    n_layers=1,
    circuit_type="Circuit_19",
)

exact = model.exact_spectrum()  # tuple of frequency arrays, one per input feature

The support is derived purely symbolically — no parameter sampling is involved. This exploits a structural property of the tree: every parameter contributes at most one sine or cosine factor per path, so the variational leaf factors are square-free monomials over \(\{1, \cos\theta_i, i\sin\theta_i\}\), which are linearly independent. A frequency is therefore absent if and only if all of its signature-grouped path weights sum to zero — an exact test, since all weights are dyadic rationals. This correctly removes frequencies whose contributions cancel identically across paths, e.g. two consecutive encodings combining into a single rotation (\(\langle Z \rangle = \cos(2x)\) has spectrum \(\{-2, 2\}\), not \(\{-2, \dots, 2\}\)).

Two methods are available:

• method="tree" (default): fully exact, but enumerates the explicit tree, whose size can grow exponentially with circuit depth.
• method="dp": merges tree nodes with identical (rotation index, observable) state, at most \(n_\text{params} \cdot 4^{n_\text{qubits}}\) states, and tracks the achievable per-feature sine/cosine counts per state as a mixed-radix bitmask. While it scales better for deep circuits and is exact per path, it cannot detect coefficients that cancel identically across paths with identical variational dependence (such as the repeated-encoding example above), where it returns a tight superset. It supports any number of input features, returning the same per-feature frequency layout as method="tree", but requires unit-magnitude input scaling: per-gate non-unit scalings such as Golomb encodings are rejected (use method="tree").

exact = model.exact_spectrum(method="dp")  # for deep circuits

Both method="dp" limitations mentioned above follow from its merged-state design, which keeps only the achievable sine/cosine counts per node and discards the per-path variational signature:

1. Cross-path cancellations are invisible. The tree method groups leaves by their variational signature and tests whether each group's weights sum to zero; the DP never stores those signatures, so it can confirm that a frequency is structurally reachable but not whether contributions from distinct paths cancel. It therefore returns a tight superset (e.g. it keeps \(\omega = 0\) for the repeated-encoding \(\langle Z \rangle = \cos(2x)\) that the tree removes).
2. Only unit-magnitude input scaling is allowed. Within a feature the DP aggregates the total sine/cosine counts \((s_f, c_f)\) over all of that feature's rotations, which is exact only when those rotations share a single angle \(x_f\). A sign flip (\(\omega = -1\), arising from Clifford commutation) is harmless because \(\cos\) is even and \(\sin\) odd, changing only the coefficient sign. A genuine per-gate scaling (as for the method="tree" case) changes the angle and would require expanding and convolving each rotation individually.

Multi-Dimensional Coefficients#

The get_spectrum method can also be used to calculate the coefficients of a model with multiple input dimensions. This feature can be enabled, by explicitly providing an encoding that supports multi-dimensional input, e.g. a list of single encodings (see Usage for details on how encodings are applied). Currently, only the FFT-based method supports this.

model = Model(
    n_qubits=4,
    n_layers=1,
    circuit_type="Circuit_19",
    random_seed=1000,
    encoding=["RX", "RY"]
)

coeffs, freqs = Coefficients.get_spectrum(model, mfs=1, mts=1
, shift=True)

psd = Coefficients.get_psd(coeffs)

Using a logarithmic color bar, one obtains the following 2D-spectrum:

Note that "X1" refers to the "RX" encoding and "X2" to the "RY" encoding.

In the multidimensional case, the freqs variable now contains the frequency indices for each dimension. This is an important detail, as due to the data_reupload argument, it is possible to have a different number of frequencies for each input dimension.

Fourier Coefficient Correlation (FCC)#

The FCC, as introduced in Fourier Fingerprints of Ansatzes in Quantum Machine Learning, is a metric that aims to predict the expected performance of an arbitrary Ansatz based on the the correlation between its Fourier modes. In this framework, the FCC for a given model can be obtained as follows:

from qml_essentials.coefficients import FCC

model = Model(
    n_qubits=6,
    n_layers=1,
    circuit_type="Hardware_Efficient",
    observables=None,
    encoding=["RY"],
)

fcc = FCC.get_fcc(
    model=model,
    n_samples=500,
    random_key=jax.random.key(1000),
)

Optionally, you can choose a different correlation method (currently "pearson", "complex_pearson", "spearman", and "covariance" are supported). The default "pearson" method preserves the existing behavior for complex coefficients by correlating stacked real and imaginary parts. The "complex_pearson" method computes a complex-valued Pearson coefficient via Hermitian normalized covariance, so the magnitude describes correlation strength and the angle describes the relative phase between coefficient vectors. The "covariance" method returns the (Hermitian, for complex coefficients) sample covariance instead of a normalized correlation; unlike the correlation methods its entries are not bounded to \([0, 1]\). When calculating the scalar FCC, calculate_fcc uses the magnitude of the fingerprint entries. Similar, other methods which require specifying n_samples (c.f. calculation of expressibility and entangling capability), methods in the FCC class take an optional parameter scale (defaults to False), which scales the number of samples depending on the number of qubits and the number of input features as \(n_\text{samples} \cdot n_\text{params} \cdot 2^{n_\text{qubits}} \cdot n_\text{features}\).

As described in our paper, the FCC is calculated as the mean of the Fourier fingerprint, which in turn can be obtained separately as follows:

fingerprint, freqs, coeffs = FCC.get_fourier_fingerprint(
    model=model,
    n_samples=500,
    random_key=jax.random.key(1000),
)

Note that actually calculating the FCC as it is shown in the paper, requires removing all the redundant entries in the fingerprint. This is implicitly done in FCC.get_fourier_fingerprint (and controlled using the trim_redundant argument), by - removing all negative frequencies (because their coefficients are complex conjugates of the positive frequencies) - removing symmetries inside the correlation matrix (the Fourier fingerprint), e.g. \(c_{0,1} = c_{1,0}\) Note that get_fcc also (by default) trims down the fingerprint before calculating the actual FCC.

When trim_redundant is enabled, the frequencies returned alongside the fingerprint are a (row_freqs, col_freqs) tuple that labels the two (trimmed) matrix axes one-to-one, so they always match the shape of the returned correlation matrix. The coefficients are returned as a matching (row_coeffs, col_coeffs) tuple whose rows align with those frequencies. The optional numerical_cap argument (passed through to get_spectrum) prunes negligible Fourier modes: coefficients below the cap are zeroed and, for a single input feature, the corresponding frequencies are dropped from the spectrum. The get_fourier_fingerprint method additionally accepts nan_to_one (default False), which replaces undefined correlation entries (NaNs, e.g. arising from constant coefficient vectors with zero variance) with one.

Both get_fcc and get_fourier_fingerprint support a weight parameter, which can be used to weight the correlation matrix, such that high-frequency components receive a lower weight. Intuitively this adresses the issue, that low frequency components have a higher impact on the mean-squared error (c.f. App. D in our paper).