Fourier Neural Operator for PDE-Based Financial Modeling

Every neural network on this blog so far has approximated a function: features in, a number out. The Fourier Neural Operator approximates an operator — a map between infinite-dimensional function spaces, where the input is a whole volatility surface and the output is a whole price surface. That is a different object, it needs different machinery, and the machinery is the point of this article: a learnable complex-valued kernel applied to the lowest Fourier modes, evaluated through an FFT round-trip in $O(n \log n)$.

Option pricing is where this lands in finance. Black-Scholes, Heston, local volatility — all PDEs, all solved today one parameter set at a time. An operator learns the entire parameter family at once, and the resulting forward pass is a single batched GPU kernel rather than a time-marching loop.

That is the promise. The honest version of this article separates the mechanics, which are solid and reproducible from the code below, from the performance claims, which in the FNO literature are reported against baselines nobody in this series would accept unexamined. The mechanics come first; the measurement agenda comes at the end, marked as unrun.

From Function Approximation to Operator Learning

Classical neural networks approximate functions: given an input $x \in \mathbb{R}^n$, they produce an output $y \in \mathbb{R}^m$. This is powerful, but fundamentally limited when the objects of interest are themselves functions. In financial PDE solving, the input is not a single number — it is a function describing initial/boundary conditions, a volatility surface, or a term structure. The output is another function: the price surface over $(S, t)$ space.

Operator learning lifts the problem to infinite-dimensional spaces. Instead of learning $f: \mathbb{R}^n \to \mathbb{R}^m$, we learn an operator:

$\mathcal{G}: \mathcal{A} \to \mathcal{U}$

where $\mathcal{A}$ and $\mathcal{U}$ are Banach spaces of functions. For option pricing, $\mathcal{A}$ might be the space of volatility surfaces $\sigma(S,t)$ and $\mathcal{U}$ the space of corresponding price surfaces $V(S,t)$.

Two architectures dominate the operator learning landscape:

1. DeepONet (Lu et al., 2021): Uses a branch network to encode the input function and a trunk network to encode the query location. The output is their inner product. Grounded in the universal approximation theorem for operators by Chen and Chen (1995).

2. Fourier Neural Operator (Li et al., 2021): Parameterizes the integral kernel in Fourier space, using the FFT for efficient global convolution. Resolution-invariant by construction.

Both are universal approximators for continuous operators, but FNO has a structural advantage for PDE problems: its spectral bias naturally captures the smooth, global structure of PDE solutions. It also has a structural disadvantage for option payoffs specifically, which we will come back to — a truncated Fourier basis and a kink at the strike are not natural allies.

FNO Architecture in Detail

The Fourier Neural Operator, introduced by Li et al. at ICLR 2021, builds on a simple but powerful observation: the Green's function (integral kernel) of many PDEs has a compact representation in Fourier space. Rather than learning a kernel in physical space — which requires $O(n^2)$ parameters for $n$ grid points — FNO learns it in frequency space with only the lowest $k_{\max}$ modes, reducing complexity to $O(n \log n)$ via the FFT.

The Iterative Architecture

An FNO consists of:

1. Lifting layer $P$: A pointwise linear map that projects the input from its original channel dimension to a higher-dimensional latent representation: $v_0(x) = P(a(x))$.

2. Fourier layers (repeated $L$ times): Each layer applies:

$v_{l+1}(x) = \sigma\left(W_l \, v_l(x) + \mathcal{K}_l(v_l)(x)\right)$

where $W_l$ is a local linear transform (pointwise $1 \times 1$ convolution) and $\mathcal{K}_l$ is a global integral operator implemented via the FFT:

$\mathcal{K}_l(v_l)(x) = \mathcal{F}^{-1}\left(R_l \cdot \mathcal{F}(v_l)\right)(x)$

Here $\mathcal{F}$ denotes the FFT, $R_l$ is a learnable complex-valued weight tensor applied to the lowest $k_{\max}$ Fourier modes, and $\sigma$ is a pointwise nonlinear activation (typically GELU).

3. Projection layer $Q$: Maps the latent representation back to the output dimension: $u(x) = Q(v_L(x))$.

Why Fourier Space?

The spectral convolution $R_l \cdot \mathcal{F}(v_l)$ is a multiplication in frequency domain, which is equivalent to a global convolution in physical space — but computed in $O(n \log n)$ instead of $O(n^2)$. This is not just an efficiency trick. PDE solutions are typically smooth and dominated by low-frequency components. By truncating to $k_{\max}$ modes, FNO acts as a learnable low-pass filter that naturally regularizes the solution and avoids high-frequency artifacts.

Crucially, FNO is claimed to be discretization-invariant: once trained on a grid of size $n$, it can be evaluated on any resolution $m \neq n$ by simply adjusting the FFT size and zero-padding or truncating the spectral weights. This zero-shot super-resolution property is unique among neural PDE solvers — and it is the first claim in this article that deserves a measurement rather than a citation. See the measurement agenda below.

The PDE We Are Learning the Operator For

The Black-Scholes PDE — derived, dissected term by term, and given with its closed-form call/put solution in The Black-Scholes Formula — is the operator target:

$\frac{\partial V}{\partial t} + \frac{1}{2}\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2} + rS\frac{\partial V}{\partial S} - rV = 0$

Everything that follows treats it as a black box with a known answer. That known answer is exactly why it is the right test case: an operator trained on finite-difference output can be scored against norm.cdf exact prices, which almost no FNO benchmark in the literature can do.

FNO Formulation

We recast the problem as operator learning. Define:

• Input function $a(S, t)$: encodes the PDE parameters. This can include the volatility surface $\sigma(S, t)$, the payoff function, and the risk-free rate as channels stacked on the $(S, t)$ grid.
• Output function $u(S, t) = V(S, t)$: the option price surface.

The FNO learns $\mathcal{G}_\theta: a \mapsto u$ from a dataset of $(a^{(i)}, u^{(i)})$ pairs generated by a traditional solver. After training, inference for any new parameter configuration is a single forward pass.

Training Data Generation

import numpy as np
from scipy.stats import norm

def black_scholes_fd(sigma, r, K, T, S_max=300, N_S=256, N_t=256):
    """Solve Black-Scholes PDE via explicit finite differences.

    NOTE: this is an interpreted double loop — the *worst* CPU baseline,
    exactly the kind called out in /en/blog/post/when-gpu-pays-off-sweep-roofline.
    It is fine for generating training data offline. It is NOT the baseline
    any speedup claim should be measured against; vectorize the inner loop
    over i (or use scipy sparse + implicit stepping) before timing anything.
    """
    dS = S_max / N_S
    dt = T / N_t
    S = np.linspace(0, S_max, N_S + 1)

    V = np.maximum(S - K, 0).astype(np.float64)

    for j in range(N_t):
        V_new = V.copy()
        for i in range(1, N_S):
            delta = (V[i+1] - V[i-1]) / (2 * dS)
            gamma = (V[i+1] - 2*V[i] + V[i-1]) / (dS**2)
            V_new[i] = V[i] + dt * (
                0.5 * sigma**2 * S[i]**2 * gamma
                + r * S[i] * delta
                - r * V[i]
            )
        V_new[0] = 0
        V_new[N_S] = S_max - K * np.exp(-r * (T - (j+1)*dt))
        V = V_new

    return S, V

For training, we sample thousands of parameter configurations — varying $\sigma \in [0.05, 0.80]$, $r \in [0.01, 0.10]$, $K \in [50, 150]$, $T \in [0.1, 2.0]$ — and solve each with the finite-difference method. The resulting dataset of input-output pairs is what the FNO learns from.

Scaling Up: The Heston Stochastic Volatility Model

Constant volatility is a known-false assumption, and why it fails — the smile, the fat tails — is the subject of the "Harsh Reality" section of The Black-Scholes Formula. Heston fixes it by making variance a second state variable:

$dS = rS\,dt + \sqrt{v}\,S\,dW_1$ $dv = \kappa(\theta - v)\,dt + \xi\sqrt{v}\,dW_2$

with $\text{Corr}(dW_1, dW_2) = \rho$. The corresponding PDE for the option price $V(S, v, t)$ is two-dimensional in space:

$\frac{\partial V}{\partial t} + \frac{1}{2}vS^2\frac{\partial^2 V}{\partial S^2} + \rho\xi vS\frac{\partial^2 V}{\partial S \partial v} + \frac{1}{2}\xi^2 v\frac{\partial^2 V}{\partial v^2} + rS\frac{\partial V}{\partial S} + \kappa(\theta - v)\frac{\partial V}{\partial v} - rV = 0$

This is where operator learning earns its keep, and the argument is structural rather than empirical. Finite-difference schemes on a $(S, v, t)$ grid scale as $O(N_S \times N_v \times N_t)$, the cross-derivative term $\partial^2 V / \partial S \partial v$ complicates the discretization, and every new parameter set $(\kappa, \theta, \xi, \rho, r)$ pays the full cost again. An operator pays it once at training time. The 2D spatial structure of the Heston PDE also maps directly onto the 2D FFT in the Fourier layers, so the architecture below generalizes with a dimension change rather than a redesign.

FNO for Option Pricing: PyTorch Implementation

Below is a complete, self-contained FNO implementation for learning the Black-Scholes operator. The architecture follows Li et al. (2021) with adaptations for the financial setting.

import torch
import torch.nn as nn
import torch.nn.functional as F
from torch.fft import rfft, irfft

class SpectralConv1d(nn.Module):
    """1D Fourier layer: spectral convolution via FFT."""

    def __init__(self, in_channels: int, out_channels: int, modes: int):
        super().__init__()
        self.in_channels = in_channels
        self.out_channels = out_channels
        self.modes = modes  # Number of Fourier modes to keep

        scale = 1.0 / (in_channels * out_channels)
        self.weights = nn.Parameter(
            scale * torch.randn(in_channels, out_channels, modes, dtype=torch.cfloat)
        )

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        batch_size = x.shape[0]

        x_ft = rfft(x, dim=-1)

        out_ft = torch.zeros(
            batch_size, self.out_channels, x_ft.size(-1),
            dtype=torch.cfloat, device=x.device
        )
        out_ft[:, :, :self.modes] = torch.einsum(
            "bix,iox->box", x_ft[:, :, :self.modes], self.weights
        )

        return irfft(out_ft, n=x.size(-1), dim=-1)


class FNOBlock(nn.Module):
    """Single Fourier Neural Operator block."""

    def __init__(self, channels: int, modes: int):
        super().__init__()
        self.spectral_conv = SpectralConv1d(channels, channels, modes)
        self.pointwise = nn.Conv1d(channels, channels, kernel_size=1)
        self.norm = nn.InstanceNorm1d(channels)

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        return F.gelu(self.norm(self.spectral_conv(x) + self.pointwise(x)))


class FNO1d(nn.Module):
    """
    Fourier Neural Operator for 1D PDE problems.

    Learns the mapping: PDE parameters -> solution function
    """

    def __init__(
        self,
        in_channels: int = 3,    # e.g., sigma(S), payoff(S), grid(S)
        out_channels: int = 1,   # V(S)
        hidden_channels: int = 64,
        modes: int = 32,
        num_layers: int = 4,
    ):
        super().__init__()
        self.lift = nn.Linear(in_channels, hidden_channels)
        self.blocks = nn.ModuleList(
            [FNOBlock(hidden_channels, modes) for _ in range(num_layers)]
        )
        self.proj = nn.Sequential(
            nn.Linear(hidden_channels, 128),
            nn.GELU(),
            nn.Linear(128, out_channels),
        )

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        x = self.lift(x)                    # -> (batch, spatial, hidden)
        x = x.permute(0, 2, 1)              # -> (batch, hidden, spatial)

        for block in self.blocks:
            x = block(x)

        x = x.permute(0, 2, 1)              # -> (batch, spatial, hidden)
        return self.proj(x)                  # -> (batch, spatial, out_channels)

Training Loop

import torch.optim as optim
from torch.utils.data import DataLoader, TensorDataset

def train_fno_black_scholes():
    device = torch.device("cuda" if torch.cuda.is_available() else "cpu")

    N_samples = 5000
    N_S = 256
    S_max = 300.0
    S_grid = np.linspace(0, S_max, N_S + 1)

    inputs, targets = [], []
    for _ in range(N_samples):
        sigma = np.random.uniform(0.05, 0.80)
        r = np.random.uniform(0.01, 0.10)
        K = np.random.uniform(50, 150)
        T = np.random.uniform(0.1, 2.0)

        _, V = black_scholes_fd(sigma, r, K, T, S_max=S_max, N_S=N_S)

        sigma_field = np.full(N_S + 1, sigma)
        payoff = np.maximum(S_grid - K, 0)
        grid_norm = S_grid / S_max

        inp = np.stack([sigma_field, payoff, grid_norm], axis=-1)
        inputs.append(inp)
        targets.append(V[:, None])

    X = torch.tensor(np.array(inputs), dtype=torch.float32)
    Y = torch.tensor(np.array(targets), dtype=torch.float32)

    dataset = TensorDataset(X, Y)
    loader = DataLoader(dataset, batch_size=64, shuffle=True)

    model = FNO1d(in_channels=3, out_channels=1, hidden_channels=64, modes=32).to(device)
    optimizer = optim.AdamW(model.parameters(), lr=1e-3, weight_decay=1e-4)
    scheduler = optim.lr_scheduler.CosineAnnealingLR(optimizer, T_max=200)

    for epoch in range(200):
        model.train()
        total_loss = 0.0
        for batch_x, batch_y in loader:
            batch_x, batch_y = batch_x.to(device), batch_y.to(device)

            pred = model(batch_x)
            loss = torch.mean(
                torch.norm(pred - batch_y, dim=1)
                / torch.norm(batch_y, dim=1).clamp(min=1e-8)
            )

            optimizer.zero_grad()
            loss.backward()
            optimizer.step()
            total_loss += loss.item()

        scheduler.step()
        if (epoch + 1) % 20 == 0:
            avg = total_loss / len(loader)
            print(f"Epoch {epoch+1:3d} | Relative L2 Loss: {avg:.6f}")

    return model

Inference: Real-Time Pricing

@torch.no_grad()
def price_option(model, sigma, K, S_grid, device="cuda"):
    """
    Price a European call for given sigma and strike.
    Returns prices for all S in S_grid — single forward pass.
    """
    S_max = S_grid[-1]
    payoff = np.maximum(S_grid - K, 0)
    grid_norm = S_grid / S_max
    sigma_field = np.full_like(S_grid, sigma)

    inp = np.stack([sigma_field, payoff, grid_norm], axis=-1)
    x = torch.tensor(inp, dtype=torch.float32).unsqueeze(0).to(device)

    pred = model(x)
    return pred.squeeze().cpu().numpy()

FNO vs. PINNs: A Practical Comparison

Physics-Informed Neural Networks and Fourier Neural Operators represent fundamentally different philosophies for solving PDEs with neural networks. Understanding their trade-offs is the practical reason to care about operator learning at all.

PINNs: Per-Instance Optimization

PINNs put the PDE residual directly in the loss — The Navier-Stokes Problem covers the mechanism, the autodiff physics_loss implementation, and the DeepMind singularity-search result. The only property that matters here is structural: a PINN is trained per parameter set. Change $\sigma$, $r$, $K$ or $T$ and you optimize again from scratch.

Strengths: No labeled data required. Directly enforces PDE structure. Works for any PDE you can write down.

Weaknesses: Must retrain for each new parameter set $(\sigma, r, K, T)$. Training involves balancing multiple loss terms (PDE residual, boundary conditions, initial conditions) which often leads to optimization pathologies. Convergence can be slow — typically 10,000–100,000 gradient steps per problem instance. Fails on multi-scale and chaotic systems where the loss landscape becomes highly non-convex.

FNO: Amortized Operator Learning

FNO learns the solution operator from data. It requires a training dataset generated by a classical solver, but once trained, it generalizes across the entire parameter space.

Strengths: Sub-millisecond inference in the published benchmarks. Generalizes to unseen parameters without retraining. Resolution-invariant — train at low resolution, evaluate at high resolution. Naturally handles smooth PDE solutions via spectral bias.

Weaknesses: Requires training data from a classical solver (chicken-and-egg problem for truly novel PDEs). Approximation error is bounded but nonzero. Less interpretable than PDE-constrained approaches. And the spectral bias that helps on smooth solutions is a liability at a payoff kink.

Head-to-Head Comparison

The table below is literature-reported, not measured here — the accuracy and speedup rows in particular come from Li et al. (2021) and follow-up work on fluid benchmarks, not from option pricing on our hardware. Read the structural rows (data needed, generalization, resolution invariance) as the reliable ones.

The Hybrid: Physics-Informed Neural Operators (PINO)

PINO (Li et al., 2024) combines both approaches. It uses the FNO architecture but augments the data-driven loss with a PDE residual term:

$\mathcal{L}_{\text{PINO}} = \mathcal{L}_{\text{data}} + \lambda \, \mathcal{L}_{\text{PDE}}$

The decomposition is the useful part. $\mathcal{L}_{\text{data}}$ anchors the operator to solver output wherever you have it; $\lambda \mathcal{L}_{\text{PDE}}$ constrains it everywhere you do not, including regions of parameter space you never sampled. For option pricing that second term is attractive for a reason specific to finance: the PDE residual is a hard constraint you can also evaluate at inference time as a self-check, which is the basis of the residual-monitoring fallback below.

Practical Considerations for Production Deployment

The generic half of this — inference latency budgets, building the data pipeline, monitoring for drift, periodic retraining — is already covered for a neural model in a trading system in the production section of DeepLOB: Deep Learning on Limit Order Books, and it applies unchanged. What follows is only what is specific to an operator.

Data Pipeline: Sampling the Parameter Space

Generating training data is the main bottleneck, and unlike a market-data model, you choose your own distribution — which means you can get it wrong in a way that is invisible in the loss. For the Black-Scholes operator with 4 parameters, 5,000–10,000 samples suffice. For Heston with 5 parameters plus a 2D spatial domain, 20,000–50,000 samples are typical. Use adaptive sampling: concentrate samples in parameter regions where the solution varies rapidly — near-the-money, short maturities, high volatility — because uniform sampling of the box spends most of its budget on deep-ITM and deep-OTM regions where the operator is nearly linear and learns fine from very few examples.

Architecture Tuning

• Modes ($k_{\max}$): Start with $N/4$ where $N$ is the spatial grid size. For $N=256$, use 32–64 modes. Too few modes lose detail near strike boundaries; too many overfit to noise.
• Layers: 4 Fourier layers are standard. Deeper networks (6–8) help for 2D problems like Heston but increase memory.
• Hidden channels: 64 for 1D Black-Scholes, 128 for 2D Heston. Scale with problem complexity.

Precision: The fp32 Question

SpectralConv1d allocates torch.cfloat — single-precision complex — and every FFT round-trip runs at that precision. This blog has already watched an fp32 financial pipeline return silent garbage in The GPU Precision Trap, where a mathematically correct prefix-sum formulation lost catastrophically at fp32 magnitudes. The analogous question here is direct: at $S \approx 100$ with prices quoted to the cent, does an fp32 spectral round-trip hold to $10^{-2}$ absolute, or does the FFT's intermediate magnitude growth eat the last significant digits?

Error Control

For production option pricing, you need error bounds you can defend.

1. Calibrated uncertainty: an ensemble standard deviation is not a coverage guarantee, and this blog has the machinery to do it properly — split conformal over a held-out parameter grid, with ACI/DtACI for the non-exchangeable case, in Conformal Prediction for Risk-Aware Position Sizing. Wrapping the FNO in split conformal over the $(\sigma, r, K, T)$ grid gives distribution-free intervals on the price.
2. Residual monitoring: Compute the PDE residual of the FNO prediction as a post-hoc check. If $\|\text{residual}\| > \epsilon$, fall back to a classical solver. This is free at inference — the residual is a finite-difference stencil on an array you already have.
3. Active learning: Route high-uncertainty inputs to the classical solver, add the results to the training set, and periodically retrain. The conformal interval width from (1) is the natural routing signal.

The Measurement Agenda

This is the part that decides whether the architecture above is worth deploying, and none of it is done yet. It is listed here rather than buried because the alternative — asserting a headline speedup — is precisely the register this blog exists to avoid.

1. Speedup is a curve, not a number. When the GPU Pays Off establishes the shape: $S(B) = aB/(O + bB)$, rising from overhead-dominated at $B=1$ to a compute-bound plateau, and it shows a headline 167x decomposing into 27x of algorithm times 6.2x of hardware. The FNO measurement must follow that template: wall-clock the forward pass at $B \in \{1, 8, 64, 512\}$ and report the whole curve. Critically, the baseline must be a vectorized, multi-core finite-difference solver — the black_scholes_fd above is an interpreted double loop, the exact "worst CPU implementation" baseline that article names, and comparing a batched GPU operator against it would produce a number that means nothing.

2. Accuracy against the closed form. This is the experiment Black-Scholes makes nearly free and that fluid-dynamics benchmarks cannot do at all: train on FD-generated data, then score against exact norm.cdf prices across the $(\sigma, r, K, T)$ box. Report the relative $L^2$, and report the error distribution — specifically near the strike and near expiry, where the solution is least smooth and the truncated spectral basis should struggle most.

3. No-arbitrage violations. A learned operator has no structural reason to respect the shape constraints a price surface must satisfy: monotonicity in $K$, convexity in $K$, and $V \geq \max(S-K, 0)$. Measuring the violation rate across the parameter box turns this article's own open question — can we guarantee no-arbitrage conditions in the learned operator? — from hand-waving into a number.

4. Discretization invariance, tested rather than asserted. Train at $N_S = 128$, evaluate at $N_S = 512$, report the error. The expected failure mode is named and specific: Gibbs ringing at the kink $\max(S-K, 0)$. A truncated Fourier basis reconstructing a function with a discontinuous first derivative oscillates around it, and zero-shot super-resolution can make that worse rather than better by exposing modes the training resolution never saw. If super-resolution degrades near the strike, that negative result is more valuable than the marketing claim it replaces — nobody in the FNO literature is pricing options with a payoff kink.

If experiments 1–4 cannot be run, this article should not ship. What remains without them is a well-written exposition of someone else's ICLR paper.

Beyond Vanilla Options

Assuming the agenda above survives contact with measurement, the FNO framework extends naturally to more complex instruments:

• American options: Add an early exercise boundary as an additional output channel. The FNO learns both the price surface and the optimal exercise boundary simultaneously.
• Barrier options: Encode barrier levels as input channels. Note that a barrier is a second discontinuity, and the Gibbs concern from the measurement agenda applies with more force, not less.
• Multi-asset baskets: Use 2D or 3D FNO for basket options on 2–3 underlyings. The curse of dimensionality is less severe than for grid-based solvers because the FNO operates on a fixed number of modes.
• Local volatility: Input the full Dupire local volatility surface $\sigma_{\text{loc}}(S,t)$ as a spatial function. This is the most natural use case for operator learning — the input is a function, not a scalar parameter.

Conclusion

The contribution of the Fourier Neural Operator is conceptual before it is computational: instead of solving one PDE at a time, you parameterize the solution operator itself, and doing that in Fourier space yields an architecture that is $O(n \log n)$, universal on continuous operators, and resolution-flexible by construction. That much follows from Li et al. (2021) and from the code above, which runs.

What does not follow is a speedup number. The literature's 100x–1000x figures are measured on fluid benchmarks against baselines whose quality is rarely stated, and this series has spent an entire article showing how much of such a headline is usually algorithm rather than hardware. The claims that would justify putting an operator in a pricing engine — cent-level fp32 accuracy, a bounded no-arbitrage violation rate, super-resolution that survives a payoff kink, and a speedup curve against a competent CPU solver — are all measurable, all cheap on Black-Scholes because the closed form exists, and all still unmeasured here.

That is the state of this draft: the mechanism is real, the promise is plausible, and the evidence is pending.

References and further reading:

• Li, Z., Kovachki, N., Azizzadenesheli, K., et al. "Fourier Neural Operator for Parametric Partial Differential Equations." ICLR 2021. arXiv:2010.08895
• Lu, L., Jin, P., Pang, G., Zhang, Z., & Karniadakis, G.E. "Learning nonlinear operators via DeepONet." Nature Machine Intelligence, 2021. doi:10.1038/s42256-021-00302-5
• Li, Z., et al. "Physics-Informed Neural Operator for Learning Partial Differential Equations." ACM/IMS Journal of Data Science, 2024. OpenReview
• neuraloperator PyTorch library: github.com/neuraloperator/neuraloperator
• Salvador, M., et al. "Neural Network Learning of Black-Scholes Equation for Option Pricing." arXiv:2405.05780
• Bai, Y., et al. "The AI Black-Scholes: Finance-Informed Neural Network." arXiv:2412.12213