Complete paper · text edition

Hyperbolic Neural Operator

Jieyuan Pei, Zhuoxuan Li, Wei Li, Haobo Zhang, Jiawei Jiang, Jianwei Zheng

All 59 pages are extracted from the author-provided PDF, including the appendix and references. Use the PDF for equation notation, figures and table layout. For an overview, see the abstract and research summary.

Hyperbolic Neural Operator
Jieyuan Pei 1 Zhuoxuan Li 2 Wei Li 1 Haobo Zhang 1 Jiawei Jiang 1 Jianwei Zheng 1 3
Abstract
Neural operators learn solution operators for para-
metric PDE families, mapping coefficients, forc-
ing fields, or geometric inputs to full solution
fields and thereby accelerating scientific computa-
tion. Transformer-based architectures offer strong
flexibility on irregular domains, but dense dot-
product attention often allocates pairwise scoring
uniformly across token pairs, neglecting that far-
field interactions in many discretized PDE ker-
nels are numerically compressible. To address
this mismatch, we draw inspiration from clas-
sical fast solvers that exploit hierarchical near–
far organization. We further observe that embed-
ding such tree-structured hierarchies in Euclidean
space incurs inherent distortion, whereas hyper-
bolic space naturally accommodates exponential
branching. Consequently, we proposeHyper-
bolic Neural Operator (HNO), which leverages
intrinsic hyperbolic geometry to instantiate a con-
tinuous Gibbs kernel based on stabilized geodesic
distances on the Lorentz hyperboloid. This design
imposes a geometric inductive bias for learnable
multi-scale near–far routing within a unified at-
tention mechanism. Empirically, HNO achieves
the lowest error among the evaluated methods on
six PDE benchmarks and two large-scale unstruc-
tured CFD tasks, reducing the mean relative ℓ2
error by up to 40% in the best evaluated setting.
Code is available in the GitHub repository.
1. Introduction
Solving partial differential equations (PDEs) is a corner-
stone of scientific and engineering computation, with broad
1College of Computer Science and Technology, Zhejiang Uni-
versity of Technology, Hangzhou, China2Tongji University, Shang-
hai, China 3Zhejiang Key Laboratory of Visual Information Intelli-
gent Processing, Zhejiang University of Technology, Hangzhou,
China. Correspondence to: Jianwei Zheng <zjw@zjut.edu.cn>.
Proceedings of the 43 rd International Conference on Machine
Learning, Seoul, South Korea. PMLR 306, 2026. Copyright 2026
by the author(s).
applications spanning electromagnetic interactions with ma-
terials and engineered structures, aerodynamic design of
vehicles and aircraft, as well as earth system and climate
modeling (Tang et al., 2025). In practice, PDEs are typically
discretized into meshes and then resolved by numerical ap-
proaches (Bonito et al., 2024); these simulations often take
several hours or even days for complex structures. Recently,
neural operators (NO) have emerged as popular surrogates
for learning solution operators of parametric PDE fami-
lies (Wu et al., 2024; Rahman et al., 2024). Given data
generated by numerical solvers, NOs map PDE inputs such
as coefficients, forcing fields, and geometric descriptors to
full solution fields, enabling orders-of-magnitude faster in-
ference (Kovachki et al., 2023), even on out-of-distribution
problems at scale (Luo et al., 2025; Hu et al., 2025).
Generally, neural operators follow the workflow of “lift-
process-project”. The input function is first lifted into a
latent representation, followed by iterative processing to
model interactions, and finally projected back to the output
function space (Kovachki et al., 2023). Currently, the schol-
arly pursuit to innovate NOs bifurcates into two branches:
spectral-based and transformer-based paradigms. Spectral
methods map inputs to representations in Fourier (Liu &
Yu, 2025), Wavelet (Hu et al., 2025), or Radon spaces (Lu
et al., 2025). These coefficient-space updates yield efficient
global coupling on regular grids; however, extending the
same mechanism to irregular geometries typically requires
additional deformation or quadrature machinery, leading to
a structural bottleneck. In contrast, transformer-based opera-
tors offer strong flexibility for irregular meshes and complex
domains. In practice, various tokenization schemes are elab-
orated, lifting functions into physics-guided slices (Luo
et al., 2025), spatio-temporal tokens (Holzschuh et al.,
2025), or geometry-aware representations in arbitrary do-
mains (Wen et al., 2025). However, a fundamental effi-
ciency gap arises from a mismatch with physical reality,
where nearby interactions require fine-grained resolution,
whereas well-separated contributions can often be summa-
rized with few coefficients (Greengard & Rokhlin, 1987).
Here and throughout,compressibilityrefers to numerical or
low-rank compressibility of far-field interactions, not fluid
incompressibility. Crucially, dense dot-product attention
overlooks this distinction, allocating pairwise scoring uni-
formly across token pairs regardless of spatial separation,
1

Hyperbolic Neural Operator
thereby rendering the computation prohibitive as the number
of tokens grows.
Classical fast solvers (Greengard & Rokhlin, 1987) exploit
this knowledge, managing near-field and far-field interac-
tions through recursive spatial partitioning, typically struc-
tured as a multilevel tree. Similarly, hierarchical matrix tech-
niques (Hackbusch, 2015) formalize this structure by com-
puting interactions within neighboring leaf clusters directly,
while approximating well-separated cluster pairs—those
satisfying an admissibility condition—via low-rank sum-
maries, as illustrated in Fig. 1B. However, embedding tree-
structured metrics in Euclidean space suffers from inherent
distortion due to a fundamental geometric mismatch: Eu-
clidean volume grows polynomially, whereas tree structures
expand exponentially (Nickel & Kiela, 2017). In this work,
we turn to hyperbolic space, which exhibits exponential vol-
ume growth relative to its radius and admits low-distortion
embeddings of tree-like metrics in low dimension (Sarkar,
2011). By contrast, Euclidean embeddings of expanding
tree metrics require substantially larger dimension or in-
cur higher distortion (Sala et al., 2018). Rather than direct
embedding of physical coordinates, we leverage intrinsic
hyperbolic geometry to modulate the interaction kernel. On
that basis, a representational bias toward near-field and far-
field distinctions naturally emerges from learned geodesic
distances. By further marrying the geometric flexibility of
transformers with the hierarchical efficiency of classical
solvers, we proposeHyperbolic Neural Operator (HNO).
Instead of treating attention as discrete token mixing, HNO
instantiates a continuous Gibbs kernel defined by negative
stabilized geodesic distances on the Lorentz hyperboloid.
HNO is inspired by the near–far organization of fast solvers,
but it does not implement multipole expansions, explicit
hierarchical trees, or error-controlled truncation; the contri-
bution is a learnable analogue of near–far routing inside the
neural operator. Our contributions are as follows.
• We propose the Hyperbolic Neural Operator (HNO),
whose core component is a stabilized Gibbs kernel
based on Lorentz geodesic distance. This design lets
scale-separation biases be learned directly within a
continuum-first framework.
• We provide a theoretical interpretation linking hyper-
bolic geometry to the near–far decomposition in nu-
merical solvers. The analysis highlights how radius–
angle coupling enables efficient routing and establishes
stability and discretization-consistency results for the
conditioned hyperbolic attention operator under stated
assumptions.
• HNO achieves the lowest error among the compared
baselines on six PDE benchmarks and two challeng-
ing unstructured-mesh CFD tasks. Empirical analysis
confirms that HNO learns to balance local precision
with global aggregation, while reducing the interaction
cost from O(N 2) to O(N M+M 2) via an M-token
interaction core.
2. Related Work
2.1. Deep Learning Solvers for PDEs
Deep learning paradigm for PDEs mainly covers instance-
specific solvers and operator learning across parametric
families. Physics-informed neural networks enforce residu-
als and boundary conditions in the continuous domain and
are commonly trained per instance (Raissi et al., 2019; Kar-
niadakis et al., 2021). Operator learning amortizes across
instances by learning maps between function spaces, includ-
ing DeepONet with branch and trunk networks (Lu et al.,
2019), kernel integral neural operators (Kovachki et al.,
2023), and transformation-based parameterizations such as
FNO (Li et al., 2021) and RNO (Lu et al., 2025).
Transformers improve long-range coupling and scalability
for operator learning. Representatively, Transolver intro-
duces Physics Attention that partitions the domain into learn-
able slices and attends to physics-aware tokens on general
geometries (Wu et al., 2024). Transolver++ further scales to
million-point geometries via highly parallel execution and a
local adaptation for industrial settings (Luo et al., 2025). Re-
cent pretrained or universal neural PDE solvers, including
OmniArch, DPOT, PDE-Transformer, and Unisolver, study
cross-family pretraining, autoregressive operator pretrain-
ing, or PDE-conditioned zero/few-shot generalization (Chen
et al., 2025; Hao et al., 2024; Holzschuh et al., 2025; Zhou
et al., 2025). These works are complementary to our setting:
HNO focuses on per-family supervised operator design and
the geometry of learned routing rather than broad cross-
family pretraining. For arbitrary domains, GAOT combines
geometry embeddings with transformer processing and mul-
tiscale encoder-decoder designs to balance accuracy and
efficiency (Wen et al., 2025). However, dense dot-product
attention in these architectures often allocates pairwise scor-
ing uniformly across token pairs, failing to exploit the com-
pressibility of far-field physical interactions.
2.2. Numerical Solvers for PDEs
Traditional numerical solvers discretize PDEs through local
schemes, including finite differences on structured grids,
finite volume methods (LeVeque, 2002), and finite element
alternatives (Brenner & Scott, 2008), with well-developed
stability and convergence theory. In nonlocal formulations
and integral-equation settings, discretization produces dense
kernel matrices. For asymptotically smooth kernels, far-field
blocks are often numerically low-rank. Fast multipole meth-
ods exploit this property to accelerate kernel summation and
2

Hyperbolic Neural Operator
Tree path
(via LCA)
A. Physical Response
(Near-far Decomposition) 
B. Hierarchical Tree
(Logical Structure) 
Global
Summary X2
X1
(space-like)
C. Hyperbolic Geometry
(Lorentz model, schematic) 
X0 (time-like)
Euclidean
crowding x
...
(space-like)
target x
Near-field
Sourse cluster C
Far-field
(summary)
ρ = ∥x − c
C ∥2
cC
ρ ≥ η-1 r
r
Well-Separated: η ≤ η0  < 1
Accuracy improves
geometrically with
expansion order  ρ
Source cluster C
Figure 1.A. Near–far decomposition.Near-field interactions require fine resolution, while far-field regions are summarized coarsely.B.
Hierarchical tree.Recursive summarization induces a multilevel cluster tree, where interactions route via lowest common ancestors
(LCA).C. Hyperbolic geometry (Lorentz model, schematic).Exponential capacity alleviates Euclidean crowding; geodesic distances
provide a proxy for hierarchy-aware routing.
matrix vector products (Greengard & Rokhlin, 1987), while
hierarchical matrix frameworks provide a unified algebraic
view and enable fast linear algebra such as approximate
factorizations and direct solutions (Hackbusch, 2015).
Kernel-independent fast multipole methods replace analytic
expansions with interpolation-based compression and re-
quire only kernel evaluations (Fong & Darve, 2009). Re-
cent work simplifies far-field compression through recur-
sive skeletonization with streamlined hierarchical data struc-
tures (Yesypenko et al., 2025), and develops translation-
invariant H2 constructions on unstructured quasi-uniform
meshes to reduce far-field storage (Börm & Henningsen,
2025). These methods are instance specific, so changes
in coefficients, geometry, or boundary conditions require
rerunning the solver and often rebuilding the compression.
3. Why Hyperbolic Geometry?
Searching for PDE solutions can be viewed as learning a
propagation operatorthat transports and aggregates infor-
mation across the domain. Classical fast solvers exploit a
near–farstructure, in which nearby interactions require fine
granularity, while contributions from well-separated regions
can be summarized with a small number of coefficients.
Fig. 1(A) visually provides this decomposition.
This compressibility stems from the fact that many PDE
kernels are asymptotically smooth (analytic away from sin-
gularities). In the far-field, contributions from a compact
source cluster C (radius r) to a distant target x (ρ≫r from
the center of C) admit truncated multipole-like expansions
with geometrically decaying truncation error; Appendix D.3
states the regularity assumptions and proof, while Green-
gard and Rokhlin (Greengard & Rokhlin, 1987) provide the
classical FMM origin.
Lemma 3.1(Far-field compressibility).An order- p far-field
expansion yields a truncation error (remainder)Rp+1(x;C)
bounded by |Rp+1(x;C)| ≤c(r/ρ) p+1, where c >0 de-
pends on the kernel and locations. This enables far-field
effects to be mediated by a small, fixed number of coeffi-
cients with controllable accuracy.
To exploit far-field compression across the whole domain,
several numerical methods apply the property recursively
across different scales. Repeated near-field resolution and
far-field summarization induce a multilevel cluster tree over
discretization elements as in Fig. 1(B). In such a hierarchy,
two locations interact strongly if they share a recent common
ancestor and weakly if their lowest common ancestor is
high in the tree. This hierarchy is induced by the chosen
clustering and discretization practices, rather than by the
physical distance of the underlying PDE domain itself.
From a geometric viewpoint, recursive near–far decompo-
sition induces a tree-structured proximity whose branching
grows exponentially with depth. Euclidean space has only
polynomial volume growth given any dimensions and suf-
fers crowding when representing such expanding structure.
The Hyperbolic space instead expands exponentially with
radius, fitting the hierarchical growth more naturally.
Lemma 3.2(Euclidean crowding and hyperbolic capacity).
Let BRm(r) and BHm(r) denote radius-r balls in Rm and
inH m with curvature−1. Their volumes satisfy
Vol

BRm(r)


∝r m,Vol

BHm(r)


∝exp

(m−1)r


,
which is consistent with representing exponentially branch-
ing structure with lower distortion in fixed dimension (Nickel
& Kiela, 2017; Ganea et al., 2018b).
As sketched in Fig. 1(C), exponential volume growth allevi-
ates Euclidean crowding and provides a geometry-aligned
proxy for hierarchy-aware routing. Taken together, far-field
compressibility and multilevel recursion motivate a scale-
tree view of operator interactions. Such scale trees sum-
marize far-field information through a bounded number of
3

Hyperbolic Neural Operator
coefficients at each level while retaining near-field fidelity.
With bounded far-field rank, dense interactions admit sub-
quadratic compression. This compression principle has long
been validated in numerical PDE solvers. However, directly
replacing dense interactions with generic linear attention can
degrade accuracy (Wu et al., 2024), since uniform lineariza-
tion does not respect near–far structure. Hyperbolic geom-
etry offers a continuous surrogate for hierarchical growth
and provides a geometry-aligned routing prior for near–far
interactions. The more formal statements and proofs are
deferred to Appendices D–G, covering cluster-tree motiva-
tion, near–far routing, Schur stability, bi-Lipschitz residual
conditions, and discretization consistency.
4. Methodology
4.1. Problem Setting
Let X=L 2(D;R da) and Y=L 2(D;R du) be Hilbert
spaces over a bounded domain D⊂R d. We focus on
learning nonlinear operators G† :X → Y that arise as
solution maps of parametric PDEs. In our PDE benchmarks,
a∈ X typically represents an input coefficient or forcing
field, andu∈ Yis the corresponding solution field.
Given ntrain observations {(aj, uj)}ntrain
j=1 with aj ∼ν
drawn i.i.d. and uj ≈ G †(aj), possibly with noise, we con-
sider a parametric family of operators
Gθ :X → Y, θ∈Θ,(1)
where Θ is a finite-dimensional parameter space and θ is
learned by minimizing the expected loss
min
θ∈Θ
Ea∼ν

ℓ

Gθ(a),G †(a)


,(2)
where ℓ is a loss functional mapping Y × Y to R. This
formulation parallels finite-dimensional supervised learn-
ing (Kovachki et al., 2023), with the key property that Θ is
shared across all discretization resolutions.
Given the problem setting, a recurring challenge ishier-
archical multiscale coupling. Recall that nearby regions
exchange fine-grained information, while distant regions in-
teract through coarser summaries. This motivates learnable
near-field and far-field routing mechanisms.
4.2. Hyperbolic Neural Operator
To instantiate the resolution-agnostic operator Gθ in Eq. (1),
we adopt the kernel-integral neural operator framework (Ko-
vachki et al., 2023). Our key departure from prior work lies
in the design of the kernel geometry. We replace the conven-
tional Euclidean dot-product similarity with a Gibbs kernel
based on stabilized hyperbolic-distance logits, which natu-
rally encodes the hierarchical near–far structure discussed
in Section 3 in a geometry-aware manner.
k (Key)
Euclidean Space
Hyperboloid Space 
    Lorentz model
time-like coordinate
Softmax Normalization
Logits Calculation
Lift:
crowding !q (Query) Minkowski
Inner Product
dot-product
Euclidean similarity
Figure 2.Hyperbolic-distance attention in the Lorentz model. We
lift query and key vectors to the hyperboloid, compute stabilized
hyperbolic distance, and convert it to attention weights with a
Gibbs kernel and softmax normalization.
Neural Operator Architecture.We follow the standard
construction comprising lifting, iterative kernel integration,
and projection.
Gθ :=Q ◦v T ◦ · · · ◦v 1 ◦ P,(3)
where P lifts input features to a higher-dimensional latent
space and Q projects back to the output space. Each iterative
update takes the form
vt+1(x) =σ t

Ltvt(x) + (Ktvt)(x) +b t(x)


,(4)
where Lt are pointwise linear operators, Kt are kernel inte-
gral operators described below, bt are learnable bias func-
tions, andσ t are activation functions.
Kernel Integral Operator.In Eq. (4), the kernel integral
operator Kt is the only term that couples distinct locations
on the continuous domain D, and thus models nonlocal
interactions. Following standard formulations (Kovachki
et al., 2023), the kernel operator acting on a vector-valued
functionv:D→R c is defined as
(Kκv)(x) =
Z
D
κ(x, y;a)v(y) dµ(y),(5)
where µ is the base measure on D, and κ(·,·;a) :D×
D→R is a scalar kernel conditioned on the input field
a. We share this scalar kernel across channels and use
learnable value and output projections WV and WO for
cross-channel interaction. Our theoretical analysis accounts
for these projections through the factor ∥WO∥op · ∥WV ∥op.
Hyperbolic Integral Operator.Standard transformer-
based neural operators parameterize κ by Euclidean dot-
product similarity, which assumes flat geometry and offers
little primitive for hierarchical near–far routing. Motivated
by Section 3, we instead instantiate κ in Eq. (5) using a
stabilized hyperbolic distance in the Lorentz model (Nickel
& Kiela, 2018). Given the input field a and the current
hidden state vt, we define query and key maps by separate
pointwise projections, written as q(x;a) :=Q θ(vt(x), a, x)
4

Hyperbolic Neural Operator
and k(y;a) :=K θ(vt(y), a, y) in Rm. In the latent cross-
attention implementation, queries are projected from latent
tokens and keys from input-point features; in latent self-
attention, queries and keys use separate latent projections.
Unlike Euclidean dot-product attention, these projections
induce adata-adaptive interaction metric, mapping (vt, a)
into a space where geodesic proximity correlates with cou-
pling strength, so the kernel provides a soft notion of near–
far interaction rather than relying on a fixed spatial partition.
The lifting to the Lorentz hyperboloid is
Π(z) :=
p
1 +∥z∥ 2, z


∈H m.(6)
The Lorentz model equips Hm with the Minkowski inner
product ⟨u, v⟩M :=−u 0v0 +Pm
i=1 uivi, which further in-
duces the geodesic distance dH(u, v) = arcosh(−⟨u, v⟩M).
For numerical stability, we score interactions by a stabilized
distance as follows.
dϵ(u, v) := arcosh

max{1 +ϵ,−⟨u, v⟩ M }


, ϵ >0.
(7)
This keeps the argument bounded away from 1; in practice,
the clamping is rarely active (Appendix E.1), so geomet-
ric properties are mainly governed by the true dH. Since
a shorter distance should yield a larger attention weight,
we form stabilized hyperbolic-distance logits as negative
distance with temperatureτ >0.
s(x, y;a) :=− 1
τ dϵ

Π(q(x;a)),Π(k(y;a))


.(8)
The resulting normalized kernel is written as follows.
κ(x, y;a) := exp(s(x, y;a))
Z(x;a) ,
Z(x;a) :=
Z
D
exp

s(x, y′;a)


dµ(y′).
(9)
Including the value and output projections, the full operator
becomesK=W O ◦ Kκ ◦W V and acts on values as
(Kv)(x) :=W O(Kκ(WV v))(x)
=W O
Z
D
κ(x, y;a)W V v(y)dµ(y). (10)
Conditioned on the query and key maps, the kernel κ(·,·;a)
is fixed and the map v7→ Kv is linear in v. See Section 4.3
for the near–far routing mechanism induced by negative
curvature.
Fast-solver view.Equations (9)–(10) formally instantiate a
conditioned integral kerneland hence a discretized kernel
matrix after Nyström quadrature, which is precisely the ob-
ject compressed by fast multipole and hierarchical-matrix
methods via near-/far-field structure. In HNO, the Lorentz
geodesic distance induces adata-adaptive near–far metric
(a) Latent Points on Hyperboloid (b) Attention Receptive Fields
LEAF (r=0.90)
MID   (r=0.82)
ROOT (r=0.19)
X0（time-like coordinate)
Figure 3.Learned multi-scale attention on Elasticity.a.Latent
tokens on the Lorentz hyperboloid.b.Corresponding receptive
fields. Tokens closer to the origin (smaller radiusr) exhibit broader
attention, while large-radius tokens attend more locally.
in feature space rather than physical coordinates, so interac-
tion strength is organized by learned geometric proximity,
yielding routing behavior analogous to near-field coupling
versus far-field summaries. This analogy is architectural and
statistical: HNO does not build an explicit FMM tree, com-
pute multipole expansions, or provide an error-controlled
truncation rule.
We approximate the integral operator with Nyström quadra-
ture on nodes {(xj, wj)}N
j=1. Let vj :=v(x j) and set
sij :=s(x i, xj;a) . The discrete weights are computed by a
weighted softmax
αij = wj exp(sij)PN
ℓ=1 wℓ exp(siℓ)
.(11)
The discrete operator then reads
(KN v)i :=W O
NX
j=1
αij WV vj.(12)
For scalability at large N, we avoid forming the dense
N×N node-to-node score matrix by introducing a com-
pression map that aggregates node features into M≪N
representative tokens. On grids and meshes, these are ob-
tained via patch tokens; on point clouds, we learn a set
of latent inducing tokens. The hyperbolic metric shapes
the non-uniform routing within this compressed interaction
pathway, paralleling hierarchical solvers that capture well-
separated interactions through low-rank surrogates rather
than exhaustive pairwise coupling.
In the analysis, we condition on the network weights and
treat ϵ >0 and τ≥τ min as fixed hyperparameters; con-
stants depend on the resulting feature-map Lipschitz bounds
Lq, Lk. Practical discretization details, including proofs of
operator bounds, are deferred to Appendix F.2.
4.3. Scale Separation and Near–Far Routing
Our attention kernel is a Gibbs kernel on hyperbolic distance
(Eq. (9)). Since the unnormalized weight is exp(−dϵ/τ),
5

Hyperbolic Neural Operator
smaller hyperbolic distances receive larger attention mass.
For 0< δ <1 , let Rδ :=τlog(1/δ) and consider the
effective interaction set {y|d ϵ(Π(q(x;a)),Π(k(y;a)))≤
Rδ}, beyond which the unnormalized mass is suppressed
by at least a factor δ. Negative curvature couples radius
and angle, yielding scale separation. That is, embeddings at
large radii interact with non-negligible weight only within
exponentially small angular apertures, whereas those near
the origin remain broadly connected. This mechanism is
visualized on Elasticity in Fig. 3.
Radius–angle coupling.We modulate attention breadth
via hyperbolic radius, with larger radii favoring locality
and smaller radii enabling broad aggregation. Let o=
(1,0, . . . ,0)be the Lorentz origin and define the radius as
r(z) :=d H

o,Π(z)


= asinh(∥z∥)
= log

∥z∥+
p
1 +∥z∥ 2

,
(13)
which grows slowly with the Euclidean norm. Hence radii
are learned implicitly through feature norms, without hard-
coding per-token scales.
The key mechanism is radius–angle coupling in hyperbolic
distance. For positive radii, write u= (u 0, u1:m) and v=
(v0, v1:m) with directions ˆu:=u 1:m/∥u1:m∥ and ˆv:=
v1:m/∥v1:m∥, and define cosθ(u, v) := ˆu⊤ˆv. For u, v∈
Hm with radiir u, rv, the hyperbolic law of cosines is
coshd H(u, v) = coshr u coshr v
−sinhr u sinhr v cosθ(u, v), (14)
so bounded distance increasingly constrains angular sep-
aration as radii grow, since the coefficient sinhr u sinhr v
scales exponentially in ru +r v. This is quantified as follows.
Proposition 4.1(Exponential angular contraction).For any
R >0 , there exists CR >0 such that if ru, rv >0 and
dH(u, v)≤R, then
θ(u, v)≤min
n
π, C R exp

− ru +r v
2
o
.(15)
Under our Gibbs kernel, the unnormalized weight decays as
exp(−dϵ/τ). Thus, through radius–angle coupling, large-
radius regions receive non-negligible attention only within
exponentially narrow cones, enabling learned separation
of near-field (local) and far-field (compressed) interactions.
Proof is deferred to Appendix E.2.
5. Experiments
We evaluate the proposed HNO on six standard PDE bench-
marks and two large-scale unstructured-mesh benchmarks,
covering a wide range of geometries and problem settings.
Benchmarks and baselines.As shown in Tab. 2, HNO
is compared with 19 baselines, covering CNN surro-
gates, frequency-domain neural operators, and transformer-
based operator learners. The suite spans point clouds,
(un)structured meshes, and regular grids under 2D and spa-
tiotemporal settings. These FNO (Li et al., 2021) and geo-
FNO (Li et al., 2023a) tasks are widely used in operator
learning. We further include two large-scale unstructured-
mesh benchmarks with about 32k nodes per sample.
Table 1.Summary of used benchmarks that cover a wide range of
geometries. #Mesh records the size of discretized meshes.
Geometry Benchmarks #Dim #Mesh
Elasticity 2D 972
Plasticity 2D+Time 3,131
Airfoil 2D 11,271
Pipe 2D 16,641
Navier-Stokes 2D+Time 4,096
Darcy 2D 7,225
ShapeNet Car 3D 32,186
AirfRANS 2D 32,000
Point
Cloud
Structured
Mesh
Regular
Grid
Unstruct.
Mesh
Implementation Setting.We follow the standard bench-
mark protocol (Wu et al., 2024) and report the mean relative
error ℓ2. For baselines, we use official implementations and
reported hyperparameters when available. HNO PDEBench
results are averaged over three runs unless otherwise noted;
baseline values follow official reports or authors’ implemen-
tations. All HNO experiments are run on a single NVIDIA
A6000 48GB GPU. The complete training details and model
configurations are provided in the Appendix I.
N = 4
N = 32
N = 128
N = 1024
Ground-Truth Hyperbolic Dot-product MLP
Figure 5.Toy multiscale tree-kernel fitting (log scale). Rows vary
tree size N; columns show the target diffusion kernel and fitted
kernels from hyperbolic, dot-product, and MLP attention.
5.1. Toy Multiscale Tree-Kernel Fitting
We construct a toy graph diffusion on a complete binary tree
of depthLwithN= 2 L leaves:
u=Ks, K ij ∝exp(−γ effdij),
NX
j=1
Kij = 1.(16)
6

Hyperbolic Neural Operator
Table 2.Main results on six standard benchmarks are reported in mean relative ℓ2 error, for which lower is better. Relative improvement
denotes the relative error reduction of our model over the second-best method, consistently computed as(e 2nd −e ours)/e2nd.
POINTCLOUD REGULARGRID STRUCTUREDMESH
OPERATOR ELASTICITY NAVIER–STOKESDARCY PLASTICITYAIRFOILPIPE
CLASSIC
UNET(2015) 0.0235 0.1982 0.0080 0.0051 0.0079 0.0065
RESNET(2016) 0.0262 0.2753 0.0587 0.0233 0.0391 0.0120
SWIN(2021) 0.0283 0.2248 0.0397 0.0170 0.0270 0.0109
DEEPONET(2019) 0.0965 0.2972 0.0588 0.0135 0.0385 0.0097
FREQUENCY
WMT (2021) 0.0359 0.1541 0.0082 0.0076 0.0075 0.0077
U-FNO (2022) 0.0239 0.2231 0.0183 0.0039 0.0269 0.0056
FNO (2021) 0.0229 0.1556 0.0108 0.0074 0.0138 0.0067
U-NO (2022) 0.0258 0.1713 0.0113 0.0034 0.0078 0.0100
F-FNO (2021) 0.0263 0.2322 0.0077 0.0047 0.0078 0.0070
LSM (2023) 0.0218 0.1535 0.0065 0.0025 0.0059 0.0050
RNO (2025) 0.0078 0.0894 0.0054 0.0012 0.0054 0.0044
TRANSFORMER
GALERKIN(2021) 0.0240 0.1401 0.0084 0.0120 0.0118 0.0098
HT-NET(2022) / 0.1847 0.0079 0.0333 0.0065 0.0059
OFORMER(2022) 0.0183 0.1705 0.0124 0.0017 0.0183 0.0168
GNOT (2023) 0.0086 0.1380 0.0105 0.0336 0.0076 0.0047
FACTFORMER(2023B) / 0.1214 0.0109 0.0312 0.0071 0.0060
ONO (2024) 0.0118 0.1195 0.0076 0.0048 0.0061 0.0052
TRANSOLVER(2024) 0.0065 0.0892 0.0057 0.0013 0.0053 0.0046
TRANSOLVER++ (2025) 0.0064 0.1010 0.0056 0.0014 0.0054 0.0042
HNO (OURS) 0.0037 0.0676 0.0045 0.0009 0.0048 0.0027
RELATIVEIMPROVEMENT +42.2% +24.2% +16.7% +25.0% +9.4% +35.7%
(a) Elasticity (b) Plasticity
 (c) Airfoil
Ground Truth Ground TruthHNO Error HNO ErrorTransolver++ Error Transolver++ Error Ground Truth HNO Error Transolver++ Error
1000
800
600
400
200
0
8
6
4
2
0
-2
-4
-6
-8
6
5
4
3
2
1
0
0.20
0.15
0.10
0.05
0.00
-0.05
-0.10
-0.15
-0.20
1.2
1.0
0.8
0.6
0.4
0.2
0.0
0.20
0.15
0.10
0.05
0.00
-0.05
-0.10
-0.15
-0.20
Figure 4.Qualitative error maps for HNO vs Transolver++.
Here s is a signal on the leaves, u is its diffused output,
and dij =d tree(i, j) is the leaf shortest-path distance deter-
mined by LCA depth. Thus K is a row-stochastic Markov
diffusion kernel whose discrete tree-distance levels induce
multiple interaction scales. Hyperbolic distance attention
preserves this nested block structure better than dot-product
and MLP controls (Fig. 5; Appendix H): the learned kernel
keeps both fine diagonal bands and coarse ancestor-level par-
titions, whereas the Euclidean controls either over-smooth
across blocks or lose the hierarchy as the tree grows. Thus
the toy study bridges the motivation and the PDE experi-
ments: it checks whether one metric can keep sharp local
bands and coherent ancestor blocks before dataset-specific
training choices enter. The scaling study in Appendix H
confirms the same pattern from small to large trees: hy-
perbolic attention consistently gives lower kernel, operator,
and row-wise KL errors than the dot-product control, while
preserving visually coherent ancestor-level blocks. This toy
problem is therefore not used as a PDE accuracy benchmark;
rather, it isolates the geometric routing bias before moving
to the full operator-learning tasks below.
5.2. Standard Benchmarks
Tab. 2 shows that HNO consistently achieves the lowest
mean relative ℓ2 errors on six standard benchmarks. For
Elasticity, nonuniform point sampling and elastic interac-
tions make multi-scale coupling difficult for Euclidean atten-
tion; for Pipe, boundary-driven flow on structured meshes
creates long-range dependencies along the channel. HNO
therefore obtains its largest gains on settings where local
resolution and global conditioning must coexist. On regu-
lar grids, HNO remains competitive, and Tab. 3 shows that
it is substantially more memory-efficient than dot–product
attention transformer baselines. Fig. 4 highlights represen-
tative error maps: HNO visibly reduces localized artifacts
compared with Transolver++, matching the numerical ad-
vantages in theRelative Improvementrow of Tab. 2. The
gains are not restricted to one discretization: HNO improves
the best competing error across point clouds, regular grids,
structured meshes, and time-dependent settings. The suite
stresses varied failure modes–irregular samples, grid fields,
boundary-driven meshes, and transient dynamics–rather
than repeating one geometry.
7

Hyperbolic Neural Operator
400
200
0
-200
-400
10
15
5
0
-5
-10
-15
HNOTransolver++RNO
(a) ShapeNet-Car: velocity error (b) ShapeNet-Car: pressure error (c-d) AirfRANS: velocity and pressure error
10
0
20
-10
-20
1.5
1.0
0.5
0.0
-0.5
-1.0
-1.5
Figure 6.Visualization of prediction errors for velocity field (left) and surface pressure (right) on ShapeNet-Car (a-b) and AirfRANS
(c-d) datasets. Rows compare HNO, Transolver++, and RNO. Blue indicates under-prediction, red indicates over-prediction, and white
represents zero error.
Figure 7.Per-sample coefficient errors: ShapeNet Car lift CL (left)
and AirfRANS dragC D (right).
5.3. Results on Large-Scale Benchmarks
We evaluate scalability on two ∼32k-node unstructured
CFD benchmarks, ShapeNet Car and AirfRANS. These set-
tings are challenging because irregular geometries require
both local boundary fidelity and long-range flow coupling.
HNO keeps the Transolver++ backbone, width/depth, to-
ken budget, radius graph, and neighbor cap, changing only
dot-product logits to stabilized hyperbolic-distance logits.
Fig. 6–7 show cleaner residual maps and lower drag/lift er-
rors than Transolver++ and RNO; details are in Appendix I.
Because graph construction, compression, and neighbor-
hood cap are fixed, the gains cannot be attributed to a larger
receptive graph or extra engineering. They instead reflect
how hyperbolic interaction geometry redistributes errors
near vehicle bodies and airfoil wakes, where boundary-layer
effects interact with global flow structures.
5.4. Ablations and Efficiency
We ablate HNO’s near–far routing on Darcy (Sec. 4.3).
Tab. 3 isolates two design choices under the same parameter
count: the interaction metric and the temperature calibra-
tion. All ablated variants use the same data split, optimizer
schedule, discretization, and backbone configuration, so the
comparison changes only the routing rule while holding
the training protocol fixed. Replacing hyperbolic distance
with Euclidean logits raises error by 10.5%, showing that
the gain comes from the geometry used to form the kernel
logits rather than from a generic Gibbs-style kernel. Fixing
Table 3.Darcy ablations and efficiency. Top: error/params; bottom:
memory/runtime.
Accuracy under matched HNO ablations
METHOD REL.ℓ 2 ↓PARAMS↓
HNO (OURS) 0.004460.0% 0.82M0.0%
−HYPERBOLIC(EUCLID) 0.00493+10.5% 0.82M+0.0%
−TEMPERATURE(FIXEDτ) 0.00464+4.0% 0.82M+0.0%
TRANSOLVER 0.00567+27.1%2.83M+246.3%
TRANSOLVER++ 0.00560+25.6%2.84M+246.3%
RNO 0.00535+19.7%2.83M+245.1%
Runtime and memory under the same benchmark protocol
METHOD
VRAM
(GB)↓
TRAIN
(H)↓
INFER
(MS/B)↓
HNO (OURS) 0.2270.0%0.730.0% 4.470.0%
−HYPERBOLIC(EUCLID) 0.227+0.0% 0.67-8.2% 4.16-6.9%
−TEMPERATURE(FIXEDτ) 0.207-8.8% 0.57-21.9% 4.49+0.4%
TRANSOLVER 2.18+860%2.76+278% 26.10+484%
TRANSOLVER++ 6.07+2574%3.12+327% 74.77+1572%
RNO 2.87+1164%5.26+621% 28.40+535%
64 80 96 120
Hidden d
Rel. 2
Hidden Dim.
2 3 4 5 6
Layers L
Layers
12 14 16 18 20
d
0.005
0.006
0.007
0.008
Q/K Dim.
Figure 8.Parameter sensitivity on Darcy for hidden dimension d,
depthL, and Q/K dimensiond H.
τ raises error by 4.0%, indicating that adaptive scale cali-
bration provides an additional but secondary benefit. The
runtime panel rules out a cost-driven explanation: HNO
stays close to its ablations in efficiency while remaining sub-
stantially lighter than transformer baselines, indicating that
accuracy is not traded for extra computation. Fig. 8 further
shows that the selected setting is not a narrow optimum, with
stable performance around moderate width/depth choices
and a compact hyperbolic Q/K dimension. Together, these
controls identify hyperbolic interaction geometry, rather
than raw capacity, extra computation, or implementation
bookkeeping, as the active ingredient.
8

Hyperbolic Neural Operator
Table 4.Matched Euclidean routing controls. Errors are in 10−3;
RF∆ is RF P90–P10 on Elasticity.
Accuracy under matched kernels(×10 −3)
Method Darcy↓Elas.↓Airf. V ol.↓Airf. Surf.↓
HNO4.455 3.688 12.009 4.674
EucDist/Dot 4.928 4.100 16.683 6.428
TempGate 4.656 5.422 71.175 8.524
RBFGate 4.768 15.020 133.168 12.560
Hierarchy induced by each scale coordinate
Method Scaleρ RF RF∆ ↑Mean RF↑
HNO radius -0.6380.275 0.331
EucDist q-norm 0.082 0.182 0.199
TempGateτ0.638 0.197 0.215
RBFGateσ0.646 0.241 0.192
5.5. Additional Analysis
Geometry Versus Euclidean Routing Controls.The
matched controls in Tab. 4 test whether HNO’s gains come
from hyperbolic geometry itself or from giving attention a
learnable notion of scale. EucDist/Dot replaces hyperbolic
distance with Euclidean dot-product-style routing, while
TempGate and RBFGate add explicit Euclidean scale vari-
ables through a learned temperature or radial-basis width.
They can still learn to make some queries broader and others
more local, as reflected by their nontrivial|ρ RF|values.
HNO nevertheless keeps the best accuracy on the two
PDEBench datasets and both AirfRANS metrics, while at-
taining the largest mean RF and RF spread on Elasticity.
This means that adding a scalar scale gate to Euclidean at-
tention is not enough to recover the same behavior. The
results support the specific role of hyperbolic radius: it acts
as a learned scale coordinate, while exponential volume
growth gives the model more room to organize multiple
global-to-local interaction ranges (Lemma D.1).
0.6
0.8
1.0H/log N
ρ = −0.61
a
-1.0
-0.5
0.0ρh(r, H)
b
1 2 3
radius r
0.2
0.4RF
ρ = −0.57
c
3 5 8 7 1 4 2 6
head
-1.0
-0.5
0.0
0.4ρh(r, RF)
d
Hyperbolic Cosine
Figure 9.Mechanism verification on Elasticity.We relate hyper-
bolic query radius r to normalized attention entropy H/logN and
physical receptive-field span RF. Larger radius indicates more
local attention, giving lower entropy and smaller RF. Panels (a,c)
show overall Spearman correlations across tokens, while panels
(b,d) show per-head correlations ρh, confirming that the radius–
locality trend holds across attention heads.
Mechanism Verification.Sec. 4.3 predicts radius-driven
scale separation under the Gibbs kernel: larger radii in-
duce local attention, while smaller radii remain more global.
On Elasticity, each query radius ri (Eq. (13)) is com-
pared with normalized attention entropy and attention span
RFi :=E j∼αi ∥xj −b i∥2. Because the diagnostic is post
hoc and the cosine baseline keeps query/key directions while
removing distance-radius coupling, it tests whether hyper-
bolic coordinates organize attention beyond ordinary angu-
lar similarity in the learned operator.
Fig. 9 shows consistent negative correlations between ra-
dius, entropy, and RF, while the cosine-similarity control
is much weaker. The head-wise trends support the same
near–far mechanism: learned radius separates broad global
aggregators from local specialists. The effect appears at the
head level rather than only after pooling, indicating multiple
operating scales rather than a single locality knob.
Table 5.Causal radius interventions on Darcy.
Method Rel.ℓ 2 ↓
HNO (baseline)0.00446
Fixed norm (Q/K) 0.03269
Shuffled norm (Q/K) 0.03328
Causal Radius Interventions.We intervene on Q/K
norms at inference while preserving angular directions and
the value pathway. Fixing the norm collapses all tokens to
a common radius; shuffling preserves the marginal norm
distribution but breaks its assignment to individual samples
and positions. Both interventions sharply degrade Darcy
accuracy (Tab. 5), supporting radius as a functional control
variable for near–far routing. Since angular directions, value
features, and under shuffling the norm histogram remain
available, the loss points to the sample- and token-specific
radial hierarchy rather than generic feature destruction.
6. Conclusion
HNO uses stabilized hyperbolic-distance logits and learned
radii to organize neural-operator attention into hierarchical
near–far interactions. This design lets the model assign dif-
ferent tokens and heads to different interaction ranges, so
local geometric details and long-range dependencies can
be handled within one attention mechanism. Across six
standard PDE benchmarks and two large-scale CFD tasks,
HNO achieves the lowest error among the evaluated meth-
ods while remaining parameter- and memory-efficient. The
toy tree-kernel study isolates the geometric bias, the Eu-
clidean routing controls show that scalar scale gates alone
do not explain the gains, and the causal-radius interventions
verify that the learned radii are functionally important. To-
gether, these results support hyperbolic organization as a
practical mechanism for neural operators on multiscale and
irregular scientific domains.
9

Hyperbolic Neural Operator
Acknowledgements
This work was supported in part by the Natural Sci-
ence Foundation of Zhejiang Province under Grant Nos.
LMS25F020006 and LZ24F030012, the China Postdoctoral
Science Foundation under Grant No. 2024M762911, and
the National Natural Science Foundation of China under
Grant Nos. 62506337 and 62276232. We thank Prof. Yan
Wang at Tsinghua AIR for providing computational support.
Impact Statement
This work develops a supervised neural-operator architec-
ture for accelerating simulations of parametric PDEs. Po-
tential benefits include faster surrogate modeling for sci-
entific and engineering design, especially on irregular or
large-scale meshes. The method does not enforce physical
constraints such as conservation laws by construction, so de-
ployment in safety-critical simulation workflows should in-
clude validation against solvers and domain-specific checks.
The experiments use scientific benchmark data and do not
involve human subjects or personal data; we do not iden-
tify societal harms beyond general risks from inaccurate
surrogate predictions if used without verification.
References
Anderson, J. W.Hyperbolic Geometry. Springer Undergrad-
uate Mathematics Series. Springer, London, 2nd edition,
2005.
Behrmann, J., Grathwohl, W., Chen, R. T. Q., Duvenaud,
D., and Jacobsen, J.-H. Invertible Residual Networks.
In Chaudhuri, K. and Salakhutdinov, R. (eds.),Proceed-
ings of the 36th International Conference on Machine
Learning, volume 97 ofProceedings of Machine Learn-
ing Research, pp. 573–582. PMLR, 09–15 Jun 2019. URL
https://proceedings.mlr.press/v97/be
hrmann19a.html.
Bonito, A., Canuto, C., Nochetto, R. H., and Veeser, A.
Adaptive finite element methods.Acta Numerica, 33:
163–485, 2024.
Börm, S. and Henningsen, J. H2-Matrices for Translation-
Invariant Kernel Functions.Engineering Analysis with
Boundary Elements, 175:106190, 2025. doi: 10.1016/j.
enganabound.2025.106190.
Brenner, S. C. and Scott, L. R.The Mathematical Theory of
Finite Element Methods. Springer, 3rd edition, 2008.
Cao, S. Choose a Transformer: Fourier or Galerkin. In
Advances in Neural Information Processing Systems, vol-
ume 34, pp. 24924–24940, 2021.
Chami, I., Ying, Z., Ré, C., and Leskovec, J. Hyperbolic
graph convolutional neural networks.Advances in neural
information processing systems, 32, 2019.
Chen, T., Zhou, H., Li, Y ., Wang, H., Gao, C., Shi, R.,
Zhang, S., and Li, J. OmniArch: Building Foundation
Model for Scientific Computing. InProceedings of the
42nd International Conference on Machine Learning, vol-
ume 267 ofProceedings of Machine Learning Research.
PMLR, 2025.
Chen, Y ., Lu, W., Xu, J., He, Y ., Li, W., and Zheng, J.
Information-Coupled Neural Operator for Computational
Mechanics and Parametric PDEs.Computer Methods in
Applied Mechanics and Engineering, 453:118851, 2026.
Fong, W. and Darve, E. The Black-Box Fast Multipole
Method.Journal of Computational Physics, 228(23):
8712–8725, 2009.
Ganea, O., Bécigneul, G., and Hofmann, T. Hyperbolic
entailment cones for learning hierarchical embeddings. In
International conference on machine learning, pp. 1646–
1655. PMLR, 2018a.
Ganea, O., Bécigneul, G., and Hofmann, T. Hyperbolic
neural networks. InAdvances in Neural Information
Processing Systems, volume 31, 2018b.
Greengard, L. and Rokhlin, V . A Fast Algorithm for Particle
Simulations.Journal of Computational Physics, 73(2):
325–348, 1987.
Gulcehre, C., Denil, M., Malinowski, M., Razavi, A., Pas-
canu, R., Hermann, K. M., Battaglia, P., Bapst, V ., Ra-
poso, D., Santoro, A., et al. Hyperbolic attention net-
works.arXiv preprint arXiv:1805.09786, 2018.
Gupta, G., Xiao, X., and Bogdan, P. Multiwavelet-Based Op-
erator Learning for Differential Equations. InAdvances
in Neural Information Processing Systems, volume 34,
pp. 24048–24062, 2021.
Hackbusch, W. Hierarchical Matrices. In Engquist, B. (ed.),
Encyclopedia of Applied and Computational Mathemat-
ics, pp. 645–647. Springer, Berlin, Heidelberg, 2015. doi:
10.1007/978-3-540-70529-1_281.
Hao, Z., Wang, Z., Su, H., Ying, C., Dong, Y ., Liu, S.,
Cheng, Z., Song, J., and Zhu, J. GNOT: A General
Neural Operator Transformer for Operator Learning. In
Proceedings of the 40th International Conference on Ma-
chine Learning, volume 202 ofProceedings of Machine
Learning Research, pp. 12556–12569. PMLR, 2023.
10

Hyperbolic Neural Operator
Hao, Z., Su, C., Liu, S., Berner, J., Ying, C., Su, H.,
Anandkumar, A., Song, J., and Zhu, J. DPOT: Auto-
Regressive Denoising Operator Transformer for Large-
Scale PDE Pre-Training. InProceedings of the 41st In-
ternational Conference on Machine Learning, volume
235 ofProceedings of Machine Learning Research, pp.
17616–17635. PMLR, 2024.
He, K., Zhang, X., Ren, S., and Sun, J. Deep Residual Learn-
ing for Image Recognition. InProceedings of the IEEE
Conference on Computer Vision and Pattern Recognition,
pp. 770–778, 2016.
He, N., Anand, R., Madhu, H., Maatouk, A., Krishnaswamy,
S., Tassiulas, L., Yang, M., and Ying, R. HELM: Hyper-
bolic Large Language Models via Mixture-of-Curvature
Experts. InAdvances in Neural Information Processing
Systems, 2025.
Holzschuh, B., Liu, Q., Kohl, G., and Thuerey, N. PDE-
Transformer: Efficient and Versatile Transformers for
Physics Simulations. InProceedings of the 42nd Interna-
tional Conference on Machine Learning, volume 267 of
Proceedings of Machine Learning Research, pp. 23562–
23602. PMLR, 2025.
Hu, P., Wang, R., Zheng, X., Zhang, T., Feng, H., Feng,
R., Wei, L., Wang, Y ., Ma, Z.-M., and Wu, T. Wavelet
Diffusion Neural Operator. InInternational Conference
on Learning Representations, 2025.
Jiang, J., Li, W., Pei, J., Zhu, J., Xu, H., Feng, Y ., and
Zheng, J. S2TA-Fuse: Semantic-Superpixel Tokenized
Attention for Spatial–Spectral Fusion.IEEE Transactions
on Geoscience and Remote Sensing, 64:1–14, 2026.
Karniadakis, G. E., Kevrekidis, I. G., Lu, L., Perdikaris,
P., Wang, S., and Yang, L. Physics-Informed Machine
Learning.Nature Reviews Physics, 3(6):422–440, 2021.
doi: 10.1038/s42254-021-00314-5.
Kovachki, N., Li, Z., Liu, B., Azizzadenesheli, K., Bhat-
tacharya, K., Stuart, A., and Anandkumar, A. Neural
Operator: Learning Maps Between Function Spaces With
Applications to PDEs.Journal of Machine Learning Re-
search, 24(89):1–97, 2023. URL https://jmlr.o
rg/papers/v24/21-1524.html.
Kress, R.Linear Integral Equations, volume 82 ofApplied
Mathematical Sciences. Springer, New York, NY , 3rd
edition, 2014. doi: 10.1007/978-1-4614-9593-2.
LeVeque, R. J.Finite Volume Methods for Hyperbolic Prob-
lems, volume 31 ofCambridge Texts in Applied Mathe-
matics. Cambridge University Press, Cambridge, 2002.
doi: 10.1017/CBO9780511791253.
Li, W., Jiang, J., Wu, J., Yu, K., and Zheng, J. LMO:
Linear Mamba Operator for MRI Reconstruction. In
Proceedings of the IEEE/CVF Conference on Computer
Vision and Pattern Recognition, pp. 5112–5122, June
2025a.
Li, W., Jiang, J., Xu, N., Cui, Y ., Li, Y ., and Zheng, J.
Spatial–Spectral Fusion Neural Operator. InIEEE Inter-
national Conference on Multimedia and Expo, 2025b.
Li, W., Zhu, J., Xu, H., Jiang, J., and Zheng, J. Spec-
Solver: Solving Spatial–Spectral Fusion via Semantic
Transformer. InProceedings of the 33rd ACM Interna-
tional Conference on Multimedia, pp. 1607–1616, 2025c.
Li, Z., Kovachki, N., Azizzadenesheli, K., Liu, B., Bhat-
tacharya, K., Stuart, A., and Anandkumar, A. Fourier
Neural Operator for Parametric Partial Differential Equa-
tions. InInternational Conference on Learning Represen-
tations, 2021.
Li, Z., Meidani, K., and Farimani, A. B. Transformer for
Partial Differential Equations’ Operator Learning.arXiv
preprint arXiv:2205.13671, 2022.
Li, Z., Huang, D. Z., Liu, B., and Anandkumar, A. Fourier
Neural Operator with Learned Deformations for PDEs
on General Geometries.Journal of Machine Learning
Research, 24(388):1–26, 2023a.
Li, Z., Shu, D., and Barati Farimani, A. Scalable Trans-
former for PDE Surrogate Modeling.Advances in Neu-
ral Information Processing Systems, 36:28010–28039,
2023b.
Liu, N. and Yu, Y . Neural Interpretable PDEs: Harmonizing
Fourier Insights with Attention for Scalable and Inter-
pretable Physics Discovery. InProceedings of the 42nd
International Conference on Machine Learning, volume
267 ofProceedings of Machine Learning Research, pp.
40085–40099. PMLR, 2025.
Liu, Q., Nickel, M., and Kiela, D. Hyperbolic Graph Neural
Networks. InAdvances in Neural Information Processing
Systems, 2019.
Liu, X., Xu, B., Cao, S., and Zhang, L. Mitigating Spectral
Bias for the Multiscale Operator Learning, 2022.
Liu, Z., Lin, Y ., Cao, Y ., Hu, H., Wei, Y ., Zhang, Z., Lin,
S., and Guo, B. Swin Transformer: Hierarchical Vision
Transformer Using Shifted Windows. InProceedings
of the IEEE/CVF international conference on computer
vision, pp. 10012–10022, 2021.
Lu, L., Jin, P., and Karniadakis, G. E. DeepONet: Learning
Nonlinear Operators for Identifying Differential Equa-
tions Based on the Universal Approximation Theorem of
Operators.arXiv preprint arXiv:1910.03193, 2019.
11

Hyperbolic Neural Operator
Lu, W., Chen, Y ., Xu, J., Li, W., Zhu, J., and Zheng, J.
Solving Partial Differential Equations via Radon Neural
Operator. InAdvances in Neural Information Processing
Systems, volume 38, 2025. URL https://papers.n
ips.cc/paper_files/paper/2025/hash/e
66233a208ef32f56df6312263239fa0-Abstr
act-Conference.html.
Luo, H., Wu, H., Zhou, H., Xing, L., Di, Y ., Wang, J., and
Long, M. Transolver++: An Accurate Neural Solver for
PDEs on Million-Scale Geometries. InProceedings of
the 42nd International Conference on Machine Learn-
ing, volume 267 ofProceedings of Machine Learning
Research, pp. 41432–41449. PMLR, 2025.
Nickel, M. and Kiela, D. Poincaré Embeddings for Learn-
ing Hierarchical Representations. InAdvances in neural
information processing systems, volume 30, 2017.
Nickel, M. and Kiela, D. Learning Continuous Hierarchies
in the Lorentz Model of Hyperbolic Geometry. InInter-
national conference on machine learning, pp. 3779–3788.
PMLR, 2018.
Rahman, A., George, R. J., Elleithy, M., Leibovici, D., Li,
Z., Bonev, B., White, C., Berner, J., Yeh, R. A., Kossaifi,
J., et al. Pretraining codomain attention neural operators
for solving multiphysics pdes.Advances in Neural Infor-
mation Processing Systems, 37:104035–104064, 2024.
Rahman, M. A., Ross, Z. E., and Azizzadenesheli, K.
U-NO: U-Shaped Neural Operators.arXiv preprint
arXiv:2204.11127, 2022.
Raissi, M., Perdikaris, P., and Karniadakis, G. E. Physics-
Informed Neural Networks: A Deep Learning Frame-
work for Solving Forward and Inverse Problems Involv-
ing Nonlinear Partial Differential Equations.Journal
of Computational Physics, 378:686–707, 2019. doi:
10.1016/j.jcp.2018.10.045.
Ronneberger, O., Fischer, P., and Brox, T. U-Net: Con-
volutional Networks for Biomedical Image Segmenta-
tion. InInternational Conference on Medical image com-
puting and computer-assisted intervention, pp. 234–241.
Springer, 2015.
Sala, F., De Sa, C., Gu, A., and Ré, C. Representation
Tradeoffs for Hyperbolic Embeddings. InInternational
conference on machine learning, pp. 4460–4469. PMLR,
2018.
Sarkar, R. Low Distortion Delaunay Embedding of Trees in
Hyperbolic Plane. InInternational symposium on graph
drawing, pp. 355–366. Springer, 2011.
Tang, Y ., Chen, R., Lou, M., Fan, J., Yu, C., Nonaka, A.,
Yao, Z., and Gao, W. Optical neural engine for solving
scientific partial differential equations.Nature Communi-
cations, 16(1):4603, 2025.
Tifrea, A., Bécigneul, G., and Ganea, O.-E. Poincaré
GloVe: Hyperbolic Word Embeddings.arXiv preprint
arXiv:1810.06546, 2018.
Tran, A., Mathews, A., Xie, L., and Ong, C. S. Fac-
torized Fourier Neural Operators.arXiv preprint
arXiv:2111.13802, 2021.
Villani, C. Stability of a 4th-order curvature condition aris-
ing in optimal transport theory.Journal of Functional
Analysis, 255(9):2683–2708, 2008.
Wen, G., Li, Z., Azizzadenesheli, K., Anandkumar, A., and
Benson, S. M. U-FNO—An enhanced Fourier neural
operator-based deep-learning model for multiphase flow.
Advances in Water Resources, 163:104180, 2022.
Wen, S., Kumbhat, A., Lingsch, L., Mousavi, S., Zhao,
Y ., Chandrashekar, P., and Mishra, S. Geometry-Aware
Operator Transformer as an Efficient and Accurate Neural
Surrogate for PDEs on Arbitrary Domains.arXiv preprint
arXiv:2505.18781, 2025.
Wu, H., Hu, T., Luo, H., Wang, J., and Long, M. Solving
High-Dimensional PDEs with Latent Spectral Models. In
Proceedings of the 40th International Conference on Ma-
chine Learning, volume 202 ofProceedings of Machine
Learning Research, pp. 37417–37438. PMLR, 2023.
Wu, H., Luo, H., Wang, H., Wang, J., and Long, M. Tran-
solver: A Fast Transformer Solver for PDEs on General
Geometries. InProceedings of the 41st International Con-
ference on Machine Learning, volume 235 ofProceed-
ings of Machine Learning Research, pp. 53681–53705.
PMLR, 2024.
Xiao, Z., Hao, Z., Lin, B., Deng, Z., and Su, H. Improved
Operator Learning by Orthogonal Attention. InProceed-
ings of the 41st International Conference on Machine
Learning, volume 235 ofProceedings of Machine Learn-
ing Research, pp. 54288–54299. PMLR, 2024.
Yang, M., Verma, H., Zhang, D. C., Liu, J., King, I., and
Ying, R. Hypformer: Exploring Efficient Transformer
Fully in Hyperbolic Space. InProceedings of the 30th
ACM SIGKDD Conference on Knowledge Discovery and
Data Mining, pp. 3770–3781, 2024. doi: 10.1145/3637
528.3672039.
Yesypenko, A., Chen, C., and Martinsson, P.-G. A Simpli-
fied Fast Multipole Method Based on Strong Recursive
Skeletonization.Journal of Computational Physics, 524:
113707, 2025.
12

Hyperbolic Neural Operator
Zheng, J., Li, W., Xu, N., Zhu, J., Lin, X., and Zhang, X.
Alias-free mamba neural operator.Advances in Neural
Information Processing Systems, 37:52962–52995, 2024.
Zhou, H., Ma, Y ., Wu, H., Wang, H., and Long, M. Uni-
solver: PDE-Conditional Transformers Towards Univer-
sal Neural PDE Solvers. InProceedings of the 42nd
International Conference on Machine Learning, volume
267 ofProceedings of Machine Learning Research, pp.
79061–79088. PMLR, 2025.
Zhu, J., Li, W., Xu, H., Jiang, J., Liu, Z., and Zheng, J.
Arbitrary-Scale Fusion Neural Operator. InProceedings
of the 33rd ACM International Conference on Multimedia,
pp. 1617–1626, 2025.
Zhu, J., Xu, H., Li, W., Jiang, J., Liu, Z., and Zheng, J.
Arbitrary-Scale Fusion Operator for High-Resolution Hy-
perspectral Imaging.IEEE Transactions on Multimedia,
2026.
13

Hyperbolic Neural Operator
Supplemental Material
Contents
Notation (Quick Reference) 16
A Hyperbolic Background 16
A.1 Hyperbolic Geometry Preliminaries (Lorentz Model and Conventions) . . . . . . . . . . . . . . . . . . . . . . . . . . 16
A.2 Hyperbolic Deep Learning Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
B Interaction Compression and Non-uniform Coupling 18
C Unsuccessful Hyperbolic Variants (What Did Not Work) 19
C.1 Family-Level Designs and Failure Analysis (Far-Field Compression View) . . . . . . . . . . . . . . . . . . . . . . . 21
D Motivation Details: Exponential Volume and Far-Field Low Rank 22
D.1 Cluster Trees in Hierarchical Solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
D.2 Exponential V olume Growth (Details for Proposition D.1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
D.3 Detailed Derivation of the Far-Field Expansion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
E Scale Separation and Near–Far Routing 26
E.1 Stabilized distanced ϵ (properties used in analysis) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
E.2 Proof of Proposition 4.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
E.3 From weight thresholds to angular cones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
F Operator Stability via Schur’s Test 27
F.1 Setup for the Hyperbolic Attention Kernel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
F.2 Boundedness via Schur’s Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
F.3 Bi-Lipschitz Residuals and Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
G Discretization Consistency 29
G.1 Setup and Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
G.2 Kernel Regularity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
G.3 Quadrature Error Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
H Toy Multiscale Tree-Kernel Fitting: Quantitative Scaling 33
I Training and Implementation Details 34
J Additional Experiments 36
J.1 Ablations and Sensitivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
J.2 Numerical Stability and Runtime Breakdown . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
J.3 Spectral Error Analysis on Darcy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
J.4 Causal Radius Interventions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
J.5 Learned Kernel Decay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
14

Hyperbolic Neural Operator
J.6 Discretization Robustness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
K Additional Visualizations 38
K.1 Qualitative Results on PDEBench . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
K.2 Direct Attention Map Visualization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
K.3 Additional Mechanism Verification Statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
15

Hyperbolic Neural Operator
Notation (Quick Reference)
Symbol Meaning
D⊂R d Spatial domain;µis the base measure onD;µ(D)its total mass.
ν,E a∼ν Data distribution over inputs and the associated expectation in the operator-learning objective.
a(·),u(·)Input field and target solution field;G † :X → Yis the ground-truth operator.
vt(·) Hidden state at layer t in the neural operator recursion; v(·) is a generic function argument of an integral
operator.
Ka(x, y),K a Abstract (data-dependent) kernel and associated integral operator in Motivation: (Kav)(x) =R
D Ka(x, y)v(y)dµ(y); later parameterized byκ(x, y;a).
α∈N d Multi-index used in the far-field/Taylor expansion (e.g.,∂ α
y and momentM α(C)).
ρ,η Far-field separation distance and separation ratio η=r/ρ (with r the source-cluster radius) in the near–far
expansion analysis.
q(·;a),k(·;a)Query/key feature maps (Euclidean, inR m) produced from inputa.
Π(·)Lorentz (hyperboloid) lift:Π(z) = (
p
1 +∥z∥ 2, z)∈H m.
⟨·,·⟩ M Minkowski inner product; induces hyperbolic distanced H.
px,m y Lifted hyperbolic query/key points:p x = Π(q(x;a)),m y = Π(k(y;a)).
r(z),r u,r v Hyperbolic radii. For z∈R m, r(z) =d H(o,Π(z)) = asinh(∥z∥) . For u, v∈H m, ru =d H(o, u) and
rv =d H(o, v), witho= (1,0, . . . ,0).
θ(u, v) For u, v∈H m with spatial parts ¯u= (u 1, . . . , um) and ¯v= (v 1, . . . , vm), define θ(u, v) =
arccos

⟨¯u,¯v⟩/(∥¯u∥ ∥¯v∥)


∈[0, π].
dϵ(·,·),τStabilized distance and temperature in logits;ϵ >0clampsarcoshnear1.
κ(x, y;a),Z(x;a)Normalized kernel and normalizer in Eq. (9). Whenais fixed, we writeκ(x, y)andZ(x)for brevity.
Kκ,KScalar kernel operator and full operator with projectionsW V , WO.
ρ,ρ h Spearman rank correlations used in mechanism analysis: ρ is computed over all heads×tokens, while ρh is
computed within each head.
{(xj, wj)}N
j=1,hNyström nodes and quadrature weights;his mesh size used in discretization bounds.
PN Sampling operator:(P N f) j =f(x j).
L,λ k Graph Laplacian (toy tree) and its eigenvalues in the discrete smoothing operator.
c0, τmin, R Lower bound on Z, temperature lower bound, and radius bound used in Schur-type and Lipschitz estimates.
A. Hyperbolic Background
A.1. Hyperbolic Geometry Preliminaries (Lorentz Model and Conventions)
Hyperbolic geometry is a standard continuous surrogate for hierarchical/tree-like structure: trees can be embedded into
hyperbolic space with low distortion, and the exponential volume growth matches branching combinatorics (Sarkar, 2011;
Nickel & Kiela, 2017; Sala et al., 2018). We summarize the Lorentz-model conventions used in this paper.
Let Rm+1 be equipped with the Minkowski bilinear form ⟨x, y⟩M :=−x 0y0 +Pm
i=1 xiyi. Define the hyperboloid model
of curvature−1:
Hm ={x∈R m+1 :⟨x, x⟩ M =−1, x 0 >0}.
For x∈H m, the tangent space is TxHm ={v∈R m+1 :⟨v, x⟩ M = 0}, and the Riemannian metric is gx(u, v) =⟨u, v⟩ M
onT xHm (positive definite).
Hyperboloid Lift.Forz∈R m, define the hyperboloid lift
Π(z) :=
p
1 +∥z∥ 2, z


∈R m+1.(17)
One verifies that ⟨Π(z),Π(z)⟩ M =−(1 +∥z∥ 2) +∥z∥ 2 =−1 and Π(z)0 =
p
1 +∥z∥ 2 >0 , so Π(z)∈H m. This
closed-form map is used to lift Euclidean feature vectors to the hyperboloid when constructing the hyperbolic attention
kernel.
Proposition A.1(Distance formula).For anyx, y∈H m, one has−⟨x, y⟩ M ≥1and
dH(x, y) = arcosh

− ⟨x, y⟩M


.
16

Hyperbolic Neural Operator
Proof. Let o= (1,0, . . . ,0)∈H m. It is standard that the proper orthochronous Lorentz group SO+(1, m) acts transitively
onH m by isometries and preserves⟨·,·⟩ M . Hence there exists an isometryg∈SO +(1, m)such thatgx=o. Then
dH(x, y) =d H(gx, gy) =d H(o, gy),⟨x, y⟩ M =⟨gx, gy⟩ M =⟨o, gy⟩ M .
By another isometry acting as a spatial rotation on (x1, . . . , xm), we may assume gy= (y ′
0, y′
1,0, . . . ,0) with y′
0 >0 and
−(y′
0)2 + (y′
1)2 =−1. Thus there exists a uniquer≥0such that(y ′
0, y′
1) = (coshr,sinhr).
Consider the curve γ: [0, r]→H m defined by γ(t) = (cosht,sinht,0, . . . ,0) . A direct computation gives
⟨γ′(t), γ′(t)⟩M = 1 , so γ is unit-speed. Moreover, γ lies in the intersection of Hm with the 2D plane spanned by o
and gy, and such intersections are geodesics in the hyperboloid model; hence γ is the unique minimizing geodesic from o to
gyand
dH(o, gy) = Length(γ) =
Z r
0
∥γ′(t)∥dt=r.
Finally,
−⟨o, gy⟩M =−⟨(1,0, . . . ,0),(y ′
0, y′
1,0, . . . ,0)⟩ M =y ′
0 = coshr= cosh(d H(o, gy)).
Applyingarcoshand using invariance undergyields the claimed formula, and also implies−⟨x, y⟩ M ≥1.
Coordinate Convention in Fig. 1C.The coordinate x0 is the time-like Minkowski coordinate imposed by ⟨x, x⟩M =−1 ;
it is not physical time. Coordinates x1, . . . , xm are space-like. Geodesics appear curved in Euclidean drawings because they
are plane sections of the hyperboloid.
A.2. Hyperbolic Deep Learning Related Work
Hyperbolic geometry has become a standard tool for representation learning when data exhibit latent hierarchies, power-law
degree distributions, or scale-free growth. Early work established that hyperbolic spaces can embed tree-like metrics with
low distortion and high parameter efficiency, popularizing Riemannian optimization in the Poincaré ball for hierarchical
embeddings (Nickel & Kiela, 2017). Subsequent work advocated the Lorentz (hyperboloid) model as a numerically stable
alternative and developed practical training recipes for learning continuous hierarchies with improved stability in high
dimensions (Nickel & Kiela, 2018). Beyond metric-preserving embeddings, hyperbolic order structures were introduced to
model asymmetric relations such as entailment and hypernymy through cone-based constructions that align with hierarchical
partial orders (Ganea et al., 2018a). These foundations motivate using hyperbolic manifolds not merely as an embedding
space, but as a geometry that matches the exponential expansion patterns of semantic and relational hierarchies.
Hyperbolic neural networks and operations.A key step from embeddings to deep learning is defining analogues of
Euclidean layers (linear maps, nonlinearities, pooling, normalization) on manifolds while controlling mapping errors
and numerical issues. Hyperbolic Neural Networks (HNN) formulated core neural operations in hyperbolic space with
Riemannian optimization, demonstrating that hyperbolic feature spaces can provide higher effective capacity than Euclidean
ones at comparable parameter budgets (Ganea et al., 2018b). In parallel, hyperbolic attention mechanisms were proposed by
replacing Euclidean dot-product similarity with hyperbolic geometry-aware compatibility functions, enabling attention to
exploit the “more room” property of hyperbolic spaces as the number of objects grows (Gulcehre et al., 2018). Hyperbolic
word embedding models further showed that distributional semantics can benefit from hyperbolic structure and can be
connected to probabilistic representations (e.g., Gaussian embeddings) while improving unsupervised performance on
similarity, analogy, and hypernymy detection (Tifrea et al., 2018). Collectively, these works establish a reusable toolbox:
(i) geometry-aware distance/similarity, (ii) manifold-compatible transformations, and (iii) training procedures that keep
optimization stable under negative curvature.
Hyperbolic graph deep learning.Graphs frequently exhibit hierarchical organization and scale-free statistics, making
them a canonical domain for hyperbolic learning. Hyperbolic Graph Convolutional Networks (HGCN) provided one of
the first inductive GCN formulations in hyperbolic space, deriving message passing and feature transformation in the
hyperboloid model and introducing trainable curvature across layers (Chami et al., 2019). Complementarily, Hyperbolic
Graph Neural Networks (HGNN) proposed architectures for learning on Riemannian manifolds via differentiable exponential
and logarithmic maps, demonstrating improved performance on benchmark graphs when hyperbolic geometry better matches
the underlying structural bias (Liu et al., 2019). These lines of work highlight a recurring theme: for relational data with
17

Hyperbolic Neural Operator
implicit hierarchy, the main gains often come from reducing representational distortion and enabling compact embeddings,
rather than merely adding depth.
Hyperbolic Transformers, scalable attention, and foundation models.More recent work pushes hyperbolic geometry
into Transformer-style architectures. Hypformer introduced a more complete hyperbolic Transformer built in the Lorentz
model, defining missing modules (e.g., linear transformations, normalization-like operations, and hyperbolic linear attention)
and demonstrating scalability on large graph benchmarks (Yang et al., 2024). At the foundation-model scale, HELM
reported billion-parameter fully hyperbolic large language models and proposed a Mixture-of-Curvature Experts mechanism
to better match heterogeneous curvature patterns in token representations, together with hyperbolic analogues of key LLM
components such as rotary positional encoding and RMS normalization (He et al., 2025). These developments indicate an
emerging direction: geometry is not only an inductive bias for hierarchy, but can also be integrated as a system-level design
choice affecting stability, scaling, and modularity in modern architectures.
B. Interaction Compression and Non-uniform Coupling
The near–far view in Section 3 treats compression as a property of the interaction operator, rather than a choice of hidden
width. In a discretization with N nodes, a standard attention layer computes a dense N×N matrix of pairwise scores
and applies a row-wise normalization to update values, resulting in O(N 2) interactions. Reducing hidden width changes
representation capacity but still retains dense coupling andO(N 2)pairwise interactions.
Interaction-level compression and complexity.HNO reduces interaction cost by routing node-to-node coupling through
M≪N representative tokens, yielding O(N M+M 2) instead of O(N 2). This complexity reduction comes from the
token bottleneck rather than the specific choice of logits. The logits mainly control non-uniform coupling and locality, while
the asymptotic cost is dictated by the N↔M routing structure. On point clouds these are learned latent tokens, while on
grids and meshes they are patch tokens. Equivalently, conditioned on routing weights, the induced node-to-node interaction
factors through anM-dimensional token space (a Nyström-like latent bottleneck).
Latent-mediated routing form.For point clouds, one layer can be written in a latent-bottleneck form. Let V∈R N×d v
denote node values and let T∈R M×d v denote latent tokens. Using row-wise normalization as in Eq. (11) with stabilized-
distance logits, we obtain two cross-attention routing matrices A∈R M×N and C∈R N×M , and a token mixing matrix
B∈R M×M . We apply
T=AV, T ′ =BT, ˜V=V+CT ′.(18)
This makes explicit that node-to-node coupling is mediated throughM tokens. Here A is row-normalized over nodes for each
token query, C is row-normalized over tokens for each node query, andB is row-normalized over tokens. When the attention
weights A, C, B are fixed, the value update is linear in V and can be written as ˜V=V+C B(AV) = (I+CBA)V .
Hence the conditioned interaction operator CBA has rank at most M. In practice, A, C, B depend on logits computed from
queries and keys, so the full layer is nonlinear and the rank statement is a conditional view of the interaction pathway. This
requires two cross-attention passes of size N×M and one token self-attention of size M×M , yielding O(N M+M 2)
interactions. For regular grids and structured meshes, patchification and unpatching are local linear maps with O(N) cost at
fixed patch size. Attention is evaluated in the patch-token space with O(M 2) cost, yielding O(N+M 2) per layer. On a 2D
grid with a fixed p×p patch size, M≈N/p 2, so the global mixing cost scales as O((N/p2)2) while patchification remains
O(N) . This refers to the global token-to-token mixing term per layer; local lifting and unlifting and pointwise updates
remainO(N). Compared to node-wiseO(N 2), this reduces cost when tokenization yieldsM≪N.
Non-uniform coupling, not uniform mixing.The bottleneck explains why the layer is cheaper, while the hyperbolic
kernel determines which interactions are emphasized. In HNO, stabilized Lorentz-distance logits define the normalized
weights in A, C, and B, inducing strongly heterogeneous routing and controllable locality across tokens and heads. This
differs from schemes whose aggregation weights are fixed or nearly uniform (e.g., constant-kernel averaging or uniform
assignments), since geometry-conditioned routing can produce both local and global receptive fields within the same
N↔M bottleneck. Uniform mixing would correspond to a row-normalized constant kernel. This yields uniform averaging
over all nodes, a row-stochastic diffusion on the node set. Empirically, attention in HNO is highly heterogeneous. For a
query token with attention weights {αj}N
j=1, define the raw entropy Hraw =− PN
j=1 αj logα j and the normalized entropy
Hnorm =H raw/logN∈[0,1] . The perplexity is P= exp(H raw) =N Hnorm, which can be interpreted as an effective
support size. Let b= PN
j=1 αjxj and define the span RF =PN
j=1 αj∥xj −b∥ 2 in physical coordinate space. On Elasticity
withN= 972points andM= 96latent queries, computed over 100 test samples,H norm ranges from0.56to1.00across
18

Hyperbolic Neural Operator
heads and latents. Equivalently, the perplexity P ranges from about 47 to 967 nodes. Moreover, hyperbolic radius controls
locality. We define the radius of a lifted embedding u∈H m as r(u) =d H(o, u) = arcosh(u0) under the Lorentz model.
The Spearman correlation between hyperbolic radius and Hnorm is about −0.61, and the Spearman correlation between
hyperbolic radius and RF is about −0.57. Tokens in the bottom 20% of radii have average perplexity about 769 and
RF≈0.38 , while tokens in the top 20% have average perplexity about 386 and RF≈0.23 . Full definitions and additional
statistics are provided in Appendix K.3. Additional visualizations are shown in Appendix K.2.
Parameter efficiency (empirical).Empirically, interaction compression is often accompanied by improved parameter
efficiency. For example on Darcy, HNO attains lower error than Transolver while using 0.82M parameters rather than
2.83M, a reduction by about3.45×.
C. Unsuccessful Hyperbolic Variants (What Did Not Work)
We additionally evaluated multiple hyperbolic operator families aimed at realizing near–far routing and far-field compression
in hyperbolic space beyond our final HNO kernel in Eq.(9). Tab. 6 reports representative Darcy runs with the same data
split and optimization setup; we additionally report executed epochs since some variants were stopped early or rerun with a
shorter budget, spanning manifold-valued backbones, learned hierarchical routers, and token-compression prototypes. These
alternatives either underperformed HNO under comparable settings or terminated early. These results delineate the design
space and are consistent with applying hyperbolic geometry to interaction weights while avoiding premature compression.
19

Hyperbolic Neural Operator
Table 6.Catalog of additional hyperbolic variants on Darcy (mean relative ℓ2 error; lower is better). We report the parameter count, the
best error observed with the same data split and optimization setup, and the executed training epochs (executed/scheduled). Status reports
the execution outcome: completed (ran the full schedule), early-stopped (terminated before schedule due to numerical or engineering
issues), underperformed (completed but remained above the anchor HNO), or crashed. Transolver and FNO are included as non-hyperbolic
references.
Variant (what changed) Params Best error Epochs Status
Anchors and references
HNO (ours)0.82M0.00446500/500 completed
Transolver(2024) 2.83M 0.00565 500/500 reference
FNO(2021) 2.38M 0.00870 500/500 reference
Irregular-to-grid patchification (sparse/soft patches)
H-SparsePatch: point →grid (bilinear), sparse patchification, patch
attention, grid→point sampling (routing logits use Euclideanℓ 2)
1.11M 0.00744 500/500 underperformed
H-SSPatch (large): soft assignment to M learned patch centers, patch
attention, decode with shared assignment weights (routing logits use
Euclideanℓ 2)
0.95M 0.01235 320/320 underperformed
Fully manifold-valued backbones
Hyperbolic Attention: manifold-valued attention backbone (Lorentz-
valued physics-attention blocks)
2.83M 0.00817 500/500 underperformed
Hyperbolic Attention (replica): second run with the same backbone 2.83M 0.00882 500/500 underperformed
HELM (native): fully manifold-valued pipeline with Lorentz lin-
ear/attention layers
2.88M 0.202 12/500 early-stopped
HELM (AMG): manifold-valued backbone with soft physics pool-
ing/unpooling
3.21M 0.0392 500/500 underperformed
Hyperbolic SSM (2D): Lorentz linear layers inside a state-space back-
bone
0.86M 0.151 44/500 early-stopped
Hyperbolic SSM (2D, 3-layer): deeper state-space backbone 1.14M 0.258 2/500 early-stopped
Scale-space / tree-routing operators with hyperbolic routing
Scale-Space Tree GAOT: hyperbolic tree router over 289 patches
(depth=4) + transformer processor
21.48M 0.0223 500/500 underperformed
Scale-Space Tree GAOT + CNN (v1): CNN spatial encoder/decoder +
contrastive regularizers
20.86M 0.0241 111/500 early-stopped
Scale-Space Tree GAOT (compact): hyperbolic encoder/decoder +
patchification + hyperbolic tree router
5.46M 0.0254 320/320 underperformed
Scale-Space Tree GAOT (compact, replica): second run with the same
configuration
5.46M 0.01049 320/320 underperformed
Scale-Space Tree GAOT (compact, lr): reduced learning rate (10−3 →
5×10 −4)
5.46M 0.01126 320/320 underperformed
Scale-Space Tree GAOT (deeper processor): increased transformer
depth (4→6layers)
8.03M 0.02154 320/320 underperformed
Scale-Space Tree GAOT (deeper processor): increased transformer
depth (4→8layers)
10.60M 0.03274 320/320 underperformed
MAGNO + hyperbolic fusion: multi-scale hyperbolic coordinate em-
bedding + neighbor aggregation + fusion
0.37M – 0/500 crashed
Other operator prototypes
Hyperbolic Slice Attention: slice-based token aggregation with
hyperbolic-distance attention
2.84M 0.00708 320/320 underperformed
H-Proto: prototype-based token compression with hyperbolic attention 1.49M 0.01047 320/320 underperformed
H-FNO: spectral backbone with hyperbolic routing/modulation 1.67M 0.00855 320/320 underperformed
20

Hyperbolic Neural Operator
C.1. Family-Level Designs and Failure Analysis (Far-Field Compression View)
Connection to Our Theory (Near–Far Compression).Section 3 frames operator learning as a near–far decomposition:
near-field interactions require high-fidelity resolution, while far-field interactions can be compressed once source/target
regions are sufficiently separated. This mirrors classical fast solvers (Greengard & Rokhlin, 1987; Hackbusch, 2015), where
far-field blocks admit low-rank/expansion-based approximations and are compressed only foradmissible(well-separated)
cluster pairs. A typical admissibility condition between two clustersCandDis
max

diam(C),diam(D)


≤ηdist(C, D), η∈(0,1),(19)
and blocks that violate (19) are refined (or treated as near-field) rather than compressed. We cite (19) only as a classical
criterion indicating when compression is provably safe in hierarchical solvers; HNO does not explicitly enforce such
a geometric test. From this view, the dominant failure mode ismisplacing compression: compression is safe only for
well-separated interactions, and compressing before near-field content is resolved introduces an irreversible bias.
A 2×2 Design Space (What We Observed).We found that the variants in Tab. 6 can be organized along two axes:
(A) how compression is applied(early N→M aggregation that averages values vs. token-mediated compression of the
interaction pathway), and(B) where hyperbolic geometry is applied(only to interactionweights/routing logits vs. to the full
value/backbonepipeline). HNO sits in the quadrant oftoken-mediated interaction compression+hyperbolic geometry only
for logits/weights. Below, we summarize each family in terms of the design axis it violates.
H-SparsePatch: Point→Grid Rasterization + Patch Attention.This variant attempted to unify irregular points and
regular grids by (i) rasterizing point features to a fixedG×G grid via bilinear scatter (producing a grid mask), (ii) patchifying
the grid with a strided convolutional patch embedding, (iii) applying patch-to-patch attention, and (iv) mapping the result
back to points via bilinear sampling. We include an Euclidean-logits ablation in this family (same architecture; routing
logits use Euclideanℓ 2) to isolate the effect of rasterization/patchification itself from the effect of geometry-aware routing.
Observed symptom:The model underperformed HNO on Darcy (Tab. 6) and on Elasticity (see additional runs below).
Why it failed:This introduces an early discretization bottleneck before any near–far routing: bilinear scatter/sampling
can alias high-frequency content and introduce boundary bias through the grid mask. Near-field details can be irreversibly
distorted, so later long-range routing cannot fully recover them; this is consistent with the systematic underperformance in
Tab. 6.
H-SSPatch: Soft Spatial Patching (Learned Centers) + Patch Attention.Instead of patchifying a grid, this variant
learned M patch centers in normalized coordinate space and produced soft assignment weights wi→m ∝exp(−∥x i −
cm∥/τ); patch tokens were aggregated as zm =P
i wi→mhi and then processed by patch-to-patch attention, followed
by decoding back to points with the same weights. As above, we report an Euclidean-logits ablation (same architecture;
routing logits use Euclidean ℓ2) to isolate the soft patching effect.Observed symptom:The method did not match HNO
on Darcy and Elasticity under comparable settings.Why it failed:This variant performs token compression at the first
step by collapsing N points into M soft patches before a near-field-preserving representation is formed. Because centers
(and a global temperature τ) impose a fixed geometric scale without an explicit admissibility test, points from distinct local
neighborhoods can be mixed across clusters in a way that is non-admissible under (19), effectively smoothing fine-scale
modes.
Fully Manifold-Valued Backbones (Hyperbolic Attention / HELM / Hyperbolic SSM (2D)).These variants moved
theentirehidden-state pipeline onto a Lorentz manifold: attention/SSM blocks required repeated exp/log maps, projection
back to the hyperboloid, and manifold-aware linear/normalization layers.Observed symptom:Some runs terminated early,
while completed runs consistently underperformed HNO even with higher parameter counts (Tab. 6).Why it failed:Our
near–far argument only requires curvature to shape interaction weights (kernels); it does not require manifold-valued value
pipelines. Making the full backbone manifold-valued adds heavy exp/log/projection overhead and radius-dependent gradient
scaling, which can destabilize optimization without improving where compression occurs.
Tree/Scale-Space Operators (Scale-Space Tree Variants).These variants replaced patchified attention with explicit
hierarchical routing: tokens were aggregated into multiscale representations and routed by learned tree/scale routers, with
hyperbolic distances used to compute routing probabilities (sometimes with additional contrastive losses).Observed
symptom:Despite being expressive (often >20 M parameters), they underperformed HNO and some configurations
21

Hyperbolic Neural Operator
terminated early (Tab. 6).Why it failed:Unlike FMM/ H-matrices—where only admissible (well-separated) blocks are
compressed (e.g., (19)) and truncation error decays with separation—these routers had no explicit structural guarantee for
what may be safely summarized. As a result, routing can summarize interactions that should remain near-field, introducing
an irreversible compression error; this is consistent with the systematic underperformance in Tab. 6. On PDEBench this is
further complicated because effective long-range coupling can be input-dependent rather than purely geometric.
Other Prototypes (Slice Assignment / Prototypes / Spectral).Hyperbolic Slice Attention used slice-based token
aggregation with hyperbolic-distance-based routing to slice centers. H-Proto compressed tokens via hyperbolic prototype
assignment and then ran attention in prototype space. H-FNO replaced the backbone with FFT-based spectral blocks and
used hyperbolic distances to modulate frequency weights (plus a hyperbolic router).Observed symptom:Hyperbolic
Slice Attention was sensitive and did not match HNO under comparable settings; H-Proto and H-FNO were stable but
substantially worse than HNO on Darcy.Why it failed:Slice/prototype assignment is an explicit token-compression step;
when assignments collapse (e.g., overly diffuse or overly peaky) or are misaligned with PDE locality, they can mix nonlocal
interactions and degrade near-field fidelity, consistent with Tab. 6. For spectral backbones, FNO already achieves global
mixing via a basis change, so adding hyperbolic routing/modulation changes interaction geometry without matching the
spectral inductive bias and typically degrades accuracy.
Additional Runs Beyond Darcy.On irregular geometries (Elasticity), we also tested (i)H-SparsePatch(grid rasterization
+ patch attention) and (ii)H-SSPatch(learned spatial patch centers + patch attention), which achieved best rel. errors
≈4.22×10 −2 and ≈1.85×10 −2 respectively (vs. HNO at 3.69×10 −3). We also prototypedH-MAGNO(MAGNO-style
neighbor aggregation + hyperbolic physics attention), but the historical run terminated early under comparable settings, so
we did not pursue it further.
Summary (Why HNO Worked While These Did Not).Overall, failures concentrated in two regions of the design space:
(i) early compression(rasterization, soft patches, slice/prototype assignment), which collapses N inputs into M summaries
before near-field structure is represented and is therefore not recoverable; and(ii) heavy or unstable routing/backbones
(explicit tree routers without an admissibility criterion such as (19), or fully manifold-valued pipelines), which add
optimization/geometry overhead without controlling whether aggregation inadvertently averages over near-field structure.
HNO keeps values in a stable Euclidean pipeline and applies hyperbolic geometry only to interactionweights, yielding
heterogeneous routing inside the token-mediated interaction core and empirically avoiding the accuracy loss observed in
early geometric-averaging variants.
D. Motivation Details: Exponential Volume and Far-Field Low Rank
D.1. Cluster Trees in Hierarchical Solvers
The multilevel hierarchy in Fig. 1Bis the same object built by fast multipole methods and hierarchical matrix methods: a
cluster treeobtained by recursively partitioning the discretization. Let {xi}N
i=1 denote discretization locations for sources
and targets.
Cluster Tree Construction in FMM.FMM builds a rooted spatial partition tree whose nodes areboxesthat contain
subsets of the points. The root is the bounding box that contains all {xi}. Each non-leaf box is subdivided into child boxes,
and points are assigned to the unique child that contains them. The recursion stops when each leaf box contains at most a
prescribed number of points. In two dimensions this subdivision is often into four child boxes and in three dimensions into
eight child boxes, but the key property is the multilevel hierarchy, not a fixed branching factor.
Index-Cluster View.Equivalently, each tree node corresponds to an index cluster C⊂ {1, . . . , N} given by the points
inside a box. The children ofCform a partition ofC, and the leaves form a disjoint partition of{1, . . . , N}.
Near–Far Separation and Routing via Ancestors.Given a cluster tree, hierarchical methods organize interactions
between two regions by scale. At fine levels, nearby clusters are handled directly to preserve near-field fidelity. At coarser
levels, interactions between well-separated clusters are approximated using a small set of coefficients, exploiting the low-rank
structure of far-field blocks for asymptotically smooth kernels. This underlies fast multipole methods and hierarchical
matrices (Greengard & Rokhlin, 1987; Hackbusch, 2015).
22

Hyperbolic Neural Operator
A useful view is that two leaf elements interact “through” their lowest common ancestor: the level of this ancestor determines
the scale at which information can be summarized without sacrificing near-field resolution. Algorithmically, this corresponds
to an upward pass that aggregates child summaries to parents, interactions between well-separated clusters at an appropriate
level, and a downward pass that propagates coarse summaries back to leaves.
D.2. Exponential Volume Growth (Details for Proposition D.1)
The following standard volume formula for hyperbolic space is a known geometric result (Anderson, 2005); we include the
short derivation for completeness.
Proposition D.1(Known hyperbolic ball volume and asymptotics).Letm≥2, letB H(r)be the geodesic ball of radiusr
inH m, and letω m−1 = vol(Sm−1). Then
vol(BH(r)) =ω m−1
Z r
0
sinhm−1(t)dt.
Moreover,
lim
r→∞
e−(m−1)r vol(BH(r)) = ωm−1
(m−1) 2 m−1 ,hencevol(B H(r))∼ ωm−1
(m−1) 2 m−1 e(m−1)r.
Proof. In geodesic polar coordinates (t, θ)∈[0,∞)×S m−1, the hyperbolic metric has the warped-product form ds2 =
dt2 + sinh2(t)dΩ 2(θ), so the Riemannian volume element is dvol = sinhm−1(t)dt dΩ . Integrating over θ and t∈[0, r]
gives the stated integral formula.
For the asymptotics, defineV(r) :=
R r
0 sinhm−1(t)dt. ThenV(r)→ ∞andV ′(r) = sinhm−1(r). By l’Hospital’s rule,
lim
r→∞
V(r)
e(m−1)r = lim
r→∞
V ′(r)
(m−1)e (m−1)r = 1
m−1 lim
r→∞
sinhr
er
m−1
.
Sincesinhr= 1
2(er −e −r), one haslim r→∞ sinhr/e r = 1/2, hence
lim
r→∞
V(r)
e(m−1)r = 1
(m−1) 2 m−1 .
Multiplying byω m−1 yields the claimed limit forvol(B H(r))and the asymptotic equivalence.
D.3. Detailed Derivation of the Far-Field Expansion
Assumption D.2(Derivative decay away fromC).Let
Bp+1(x;C) := sup
z∈co(C∪{c})
max
|α|=p+1
∂α
y Ka(x, z)
.
For each input fielda, there existsA p+1(x;a)∈(0,∞)such that for every bounded measurableCand everyc∈R d,
Bp+1(x;C)≤ Ap+1(x;a)
ρp+1 whereρ= dist

x, co(C∪ {c})


.
Lemma D.3(Taylor far-field moment expansion).Assume that for fixed x, the map y7→K a(x, y) is (p+ 1) -times
continuously differentiable on a neighborhood ofco(C∪ {c}) and satisfies the derivative decay condition in Assumption D.2.
When the source region C is well-separated from the evaluation point x and the separation ratio satisfies η=r/ρ <1 , the
kernel integral admits a multipole expansion:
Z
C
Ka(x, y)s(y)dy=
X
|α|≤p
∂α
y Ka(x, c)
α! Mα(C) +R p+1(x;C), (20)
where Mα(C) =
R
C(y−c) αs(y)dy are aggregated moments, c is the cluster center, r is the cluster radius, and ρ=
dist

x, co(C∪ {c})


. Moreover, the remainder satisfies|Rp+1(x;C)| ≤ C(x, p;a)η p+1, where C(x, p;a) is derived below.
23

Hyperbolic Neural Operator
Setup on Measurability, Norms, and Separation.Let C⊂R d be Lebesgue measurable and bounded, let c∈R d be
fixed, and let s∈L 1(C). Fix a norm ∥ · ∥ on Rd. Let Ceq ≥1 be such that ∥v∥∞ ≤C eq∥v∥ for all v∈R d. Such a constant
exists by norm equivalence. We define
r= sup
y∈C
∥y−c∥<∞, ρ= dist

x, co(C∪ {c})


>0, η=r/ρ <1.
If c∈C , which holds for typical choices of region centers, then r≤diam(C) and hence η≤diam(C)/ρ . More generally
r≤diam(C) + dist(c, C).
Assume that for fixed x, the map y7→K a(x, y) is (p+ 1) -times continuously differentiable on an open convex set U⊂R d
such that
co(C∪ {c})⊂U.
In particular, the segment[c, y] ={c+t(y−c) :t∈[0,1]} lies in U for each y∈C . This ensures that multivariate Taylor’s
theorem with integral remainder along[c, y]applies and that the derivative envelopeB p+1(x;C)defined below is finite.
Taylor Expansion with Integral Remainder.For a multi-index α∈N d, write |α|= P
i αi, α! =Q
i αi!, (y−c) α =Q
i(yi −c i)αi, and∂ α
y =Q
i ∂αi
yi . For eachy∈C, Taylor’s theorem aroundcyields
Ka(x, y) =
X
|α|≤p
1
α! ∂α
y Ka(x, c) (y−c) α +
X
|α|=p+1
p+ 1
α! (y−c) α
Z 1
0
(1−t) p ∂α
y Ka

x, c+t(y−c)


dt.(21)
Define the moments (well-defined sinces∈L 1(C)andCis bounded):
Mα(C) =
Z
C
(y−c) α s(y)dy,|α| ≤p.
Moment Expansion Justified by Tonelli.Let
I(x;C) =
Z
C
Ka(x, y)s(y)dy.
Multiplying (21) bys(y)and integrating overC, the finite sum over|α| ≤pis immediate.
For later use, we define the order(p+ 1)derivative envelope over the relevant compact set.
Bp+1(x;C) = sup
z∈co(C∪{c})
max
|α|=p+1
∂α
y Ka(x, z)
 <∞.
Moreover, for|α|=p+ 1and(y, t)∈C×[0,1],
|(y−c) α| ≤ ∥y−c∥ p+1
∞ ≤C p+1
eq ∥y−c∥ p+1,
∂α
y Ka

x, c+t(y−c)

 ≤B p+1(x;C),
so the integrand is dominated by C p+1
eq Bp+1(x;C)∥y−c∥ p+1|s(y)|(1−t) p, which is integrable since C is bounded,
s∈L 1(C), and
R 1
0 (1−t) pdt <∞. Hence Tonelli’s theorem applies to the remainder term. Therefore
I(x;C) =
X
|α|≤p
1
α! ∂α
y Ka(x, c)M α(C) +R p+1(x;C),(22)
where the truncation remainder admits the exact representation
Rp+1(x;C) =
X
|α|=p+1
p+ 1
α!
Z
C
(y−c) αs(y)
Z 1
0
(1−t) p ∂α
y Ka

x, c+t(y−c)


dt dy.(23)
24

Hyperbolic Neural Operator
Deterministic Remainder Bound.Using (23), the bound |(y−c) α| ≤ ∥y−c∥ |α|
∞ ≤C |α|
eq ∥y−c∥ |α|, and
R 1
0 (1−t) pdt=
1
p+1,
|Rp+1(x;C)| ≤
X
|α|=p+1
p+ 1
α!
Z
C
|(y−c) α| |s(y)|
Z 1
0
(1−t) p∂α
y Ka(x, c+t(y−c))
 dt dy
≤
X
|α|=p+1
p+ 1
α! Bp+1(x;C)C p+1
eq
Z
C
∥y−c∥ p+1|s(y)|
Z 1
0
(1−t) p dt dy
=B p+1(x;C)C p+1
eq
 X
|α|=p+1
1
α!
Z
C
∥y−c∥ p+1|s(y)|dy
≤B p+1(x;C)C p+1
eq
 X
|α|=p+1
1
α!

rp+1
Z
C
|s(y)|dy.(24)
The combinatorial factor can be computedexactly:
X
|α|=n
1
α! = dn
n! , n∈N,
since the coefficient of tn in exp(t)d = Qd
i=1 exp(t) equalsP
|α|=n 1/α! and also equals dn/n!. Applying this with
n=p+ 1to (24) gives the explicit bound
|Rp+1(x;C)| ≤ dp+1
(p+ 1)! Bp+1(x;C)C p+1
eq rp+1 ∥s∥L1(C).(25)
Separation-Rate Form (ηp+1) Under a Kernel Regularity Assumption.The step from (25) to an ηp+1 rate requires the
distance-to-singularity derivative control stated in Assumption D.2.
Theorem D.4(Far-field remainder bound).Under the setup of Lemma D.3, with C bounded, s∈L 1(C), ρ >0 ,
η=r/ρ <1, and Assumption D.2, the truncation remainder in(22)satisfies
|Rp+1(x;C)| ≤ C(x, p;a)η p+1,(26)
whereC(x, p;a) := dp+1
(p+1)! Ap+1(x;a)C p+1
eq ∥s∥L1(C) andη:=r/ρ.
Proof.Combining (25) withr=ηρand Assumption D.2 yields the result.
Corollary D.5(Far-field low rank: rank–accuracy tradeoff).Let rp :=|{α∈N d :|α| ≤p}|=
p+d
d


. Define the truncated
far-field operator (cf. Eq.(20))
(Tps)(x) :=
X
|α|≤p
∂α
y Ka(x, c)
α! Mα(C), M α(C) =
Z
C
(y−c) αs(y)dy.
Then Tp factors through the rp moments {Mα(C)} |α|≤p and is therefore a rank-rp separated approximation of the far-field
block. Moreover, the approximation error is controlled by Theorem D.4:

Z
C
Ka(x, y)s(y)dy−(T ps)(x)
 ≤ C(x, p;a)η p+1.
In particular, for any finite set of target points{xi}, the induced interaction matrix between {xi} and sources in C admits a
rank-rp approximation with accuracy governed byη p+1.
Proof. The representation is a sum of rp separable terms indexed by |α| ≤p , so it factors through an rp-dimensional
moment vector. The error bound is exactly (26).
25

Hyperbolic Neural Operator
Table 7.Empirical clamp activation in Eq. (7) on benchmark evaluations using our best checkpoints. We report the count of pairwise
distance evaluations where−⟨u, v⟩M <1+ϵ (i.e., the clamping path is taken), along with the minimum observed marginmin(−⟨u, v⟩M −
1).
Datasetϵ#pairs Clamp countmin(−⟨u, v⟩ M −1)
Elasticity10 −4 2.38×10 8 0 0.370
Navier–Stokes10 −6 2.10×10 9 0 0.717
Darcy10 −6 2.67×10 8 0 0.511
Plasticity10 −4 8.95×10 9 0 0.183
Airfoil10 −6 1.52×10 10 0 0.186
Pipe10 −6 1.09×10 10 0 0.567
RemarkD.6 (Constant dependencies).The constant C(x, p;a) in Theorem D.4 depends only on d, p, the derivative-decay
envelope Ap+1(x;a) , the norm-equivalence constant Ceq, and ∥s∥L1(C). No dependence on the shape of C enters beyond η
and∥s∥ L1(C).
RemarkD.7 (When Assumption D.2 holds).For classical PDE kernels (Laplace/Helmholtz/heat) away from their singulari-
ties, bounds of the formmax|α|=p+1 |∂α
y K(x, y)| ≤C p+1∥x−y∥ −(p+1) hold on {∥x−y∥ ≥ρ} , implying Assumption D.2
with ρ= dist

x, co(C∪ {c})


. For general data-dependent kernels Ka, Assumption D.2 should be treated as an explicit
regularity condition. This assumption is used only to formalize the far-field motivation from classical fast solvers; our
subsequent stability (Appendices E and F.2) and discretization (Appendix G) analyses do not rely on it.
E. Scale Separation and Near–Far Routing
E.1. Stabilized distanced ϵ (properties used in analysis)
RemarkE.1 ( dϵ is not a metric).The clamping in Eq. (7) can break the triangle inequality, so dϵ is generally not a metric on
Hm. In this paper, we only use dϵ as a stabilized similarity scale for logits, together with the monotonicity and Lipschitz
properties below. Empirically, the clamp is never activated on our benchmark evaluations (Tab. 7), so the behavior is
governed by the true hyperbolic distance dH. In other words, the clamp is a training-time numerical safeguard (to avoid the
singular derivative ofarcoshat1), and it is inactive at inference for our best checkpoints.
Lemma E.2(Monotonicity of dϵ and logits).Let g(t) = arcosh(max{t,1 +ϵ}) with ϵ >0 . Then g is nondecreasing,
hence dϵ(u, v) =g(−⟨u, v⟩ M) is nondecreasing in −⟨u, v⟩M . Consequently, for any fixed τ >0 , the logit s=−d ϵ/τ is
nonincreasing in−⟨u, v⟩ M .
Proof. The map t7→max{t,1 +ϵ} is nondecreasing, and arcosh is nondecreasing on [1,∞) , so their composition g is
nondecreasing. The remaining statements follow by composition with−⟨u, v⟩ M .
Lemma E.3(Global Lipschitz bound induced by clamping).Let g(t) = arcosh(max{t,1 +ϵ}) with ϵ >0 . Then g is
globally Lipschitz onRwith
Lip(g)≤ 1p
(1 +ϵ) 2 −1
= 1p
ϵ(2 +ϵ)
.
In particular, for anyu 1, u2, v∈H m,
dϵ(u1, v)−d ϵ(u2, v)
 ≤ 1p
ϵ(2 +ϵ)
⟨u1 −u 2, v⟩M
 ≤ ∥v∥2p
ϵ(2 +ϵ)
∥u1 −u 2∥2,
and similarly in the second argument. On bounded subsets of Hm, this yields a uniform Euclidean Lipschitz constant and is
the only place where the factor 1/
p
ϵ(2 +ϵ) enters; Lemma G.1 translates this into Lipschitzness of the kernel in physical
coordinates.
Proof. For t >1 , d
dt arcosh(t) = 1√
t2−1. Since max{t,1 +ϵ} ≥1 +ϵ , the derivative is bounded by 1/
p
(1 +ϵ) 2 −1
wherever it exists, and g is constant on (−∞,1 +ϵ] . Hence g is globally Lipschitz with the stated constant. The final
inequalities use|⟨a, b⟩ M | ≤ ∥a∥ 2∥b∥2.
26

Hyperbolic Neural Operator
E.2. Proof of Proposition 4.1
On the Lorentz hyperboloid, write u= (u 0,¯u)and v= (v 0,¯v). For u∈H m, we have u0 = coshr u and ∥¯u∥= sinhru,
wherer u =d H(o, u), and similarly forv. Assumer u, rv >0and letθ=θ(u, v)satisfycosθ=⟨¯u,¯v⟩/(∥¯u∥∥¯v∥).
The hyperbolic law of cosines gives
coshd H(u, v) = coshr u coshr v −sinhr u sinhr v cosθ= cosh(r u −r v) + 2 sinhru sinhr v sin2 θ
2 .
Assumed H(u, v)≤R. Then
2 sinhru sinhr v sin2 θ
2 ≤coshR−cosh(r u −r v)≤coshR.
Forr≥1,sinhr≥ 1
2(er −1)≥ 1
2 er−1, hence forr u, rv ≥1,
sinhr u sinhr v ≥ 1
4 eru+rv−2.
Therefore,
sin2 θ
2 ≤ coshR
2 sinhru sinhr v
≤2e 2 coshR·e −(ru+rv).
Usingsint≥2t/πfort∈[0, π/2](concavity ofsin) yieldssin(θ/2)≥θ/πforθ∈[0, π], hence
θ≤π
√
2e2 coshRexp

− ru +r v
2

=:C R exp

− ru +r v
2

,
which proves the claim for ru, rv ≥1 . If instead min{ru, rv}<1 , then |ru −r v| ≤d H(u, v)≤R implies ru +r v ≤R+2 ,
so exp(−(ru +r v)/2)≥exp(−(R+ 2)/2) . Since θ∈[0, π] and CR exp(−(R+ 2)/2)≥π (using coshR≥e R/2), the
same bound holds for all radii.
E.3. From weight thresholds to angular cones
Corollary E.4(Thresholded attention concentrates in a hyperbolic ball and an angular cone).Fix a query point u∈H m, a
measurable key map y7→v(y)∈H m, and τ >0 . Let dmin := inf y∈D dϵ(u, v(y)) and define the (relative) δ-threshold set
forδ∈(0,1),
Nδ(u) :=
n
y∈D: exp

−d ϵ(u, v(y))/τ


≥δexp

−d min/τ

o
.
Equivalently, if κu(y)∝exp(−d ϵ(u, v(y))/τ) denotes the normalized attention weight for fixed u, then Nδ(u) ={y:
κu(y)≥δsup y′ κu(y′)}. Then Nδ(u)⊆ {y∈D:d ϵ(u, v(y))≤d min +τlog(1/δ)} . Moreover, since dH(u, v)≤
dϵ(u, v), one also has
Nδ(u)⊆
n
y∈D:d H(u, v(y))≤d min +τlog(1/δ)
o
,
and Proposition 4.1 implies that points in Nδ(u) lie in an angular cone whose aperture decays as exp(−(ru +r v(y))/2),
up to a constant depending ond min +τlog(1/δ).
Proof. The first inclusion is immediate: exp(−dϵ/τ)≥δexp(−d min/τ) is equivalent to dϵ ≤d min +τlog(1/δ) . The
second inclusion uses dH ≤d ϵ, which follows from max{1 +ϵ,−⟨u, v⟩ M } ≥ −⟨u, v⟩ M and monotonicity of arcosh.
Applying Proposition 4.1 withR=d min +τlog(1/δ)yields the cone bound.
F. Operator Stability via Schur’s Test
F.1. Setup for the Hyperbolic Attention Kernel
LetD⊂R d be a bounded domain with measureµ. We distinguish two operator levels.
Scalar Kernel Operator.The scalar attention operatorK κ :L 2(D)→L 2(D)is defined by
(Kκv)(x) =
Z
D
κ(x, y;a)v(y)dµ(y),(27)
whereκ(x, y;a)is the normalized hyperbolic attention kernel, a scalar.
27

Hyperbolic Neural Operator
Full Operator with Projections.With learnable projectionsW O, WV , the full operator is
K:=W O ◦ Kκ ◦W V .
All Schur bounds are first derived forK κ, then∥K∥ ≤ ∥W O∥op∥WV ∥op∥Kκ∥.
Vector-Valued and Multi-Head Extension.The analysis in the main text focuses on a scalar kernel operator Kκ and
accounts for value/output projections by operator norms. The same reasoning extends to standard multi-head attention by
viewing heads as a direct sum of bounded operators.
Lemma F.1(Multi-head operator norm bound).Let H≥1 . For each head h∈ {1, . . . , H} , let K(h)
κ :L 2(D)→L 2(D)
be a scalar-kernel operator with ∥K(h)
κ ∥L2→L2 ≤C h, and denote by the same symbol its componentwise extension to
vector-valued functions. Let W (h)
V be the value projection for head h and let WO be the (pointwise) output projection
applied after concatenation. Define the multi-head layer on vector-valued functions by
KMH(v) :=W O
h
K(1)
κ (W (1)
V v), . . . ,K(H)
κ (W (H)
V v)
i
,
where[·, . . . ,·]denotes concatenation over heads. Then
∥KMH∥L2(D;Rdv )→L2(D;Rdv ) ≤ ∥W O∥op
 HX
h=1

Ch∥W (h)
V ∥op

21/2
≤ ∥W O∥op
√
Hmax
h

Ch∥W (h)
V ∥op


.
Proof.By definition and Cauchy–Schwarz in the direct-sum Hilbert space,
HX
h=1
∥K(h)
κ (W (h)
V v)∥2
L2 ≤
HX
h=1
C2
h∥W (h)
V ∥2
op ∥v∥2
L2 .
Concatenation preserves the squared norm by summation over heads, and applying WO pointwise contributes a factor
∥WO∥op.
F.2. Boundedness via Schur’s Test
We show boundedness ofK κ for the continuum normalized kernel.
F.2.1. CONTINUUMNORMALIZEDKERNEL
For the continuum normalized kernel (Eq. 9):
κ(x, y) = exp(−dϵ(Π(q(x)),Π(k(y)))/τ)
Z(x) , Z(x) :=
Z
D
exp(−dϵ(Π(q(x)),Π(k(y ′)))/τ)dµ(y ′).
Proposition F.2(Continuum Boundedness).The continuum normalized kernel satisfies:
1.κ(x, y)≥0and
R
D κ(x, y)dµ(y) = 1for allx(normalization).
2. IfZ(x)≥c 0 >0uniformly, then
R
D κ(x, y)dµ(x)≤µ(D)/c 0.
By Schur’s test (Kress, 2014) withC1 = 1andC 2 =µ(D)/c 0:
∥Kκ∥L2→L2 ≤
p
µ(D)/c0,∥K∥ L2→L2 ≤ ∥W O∥op∥WV ∥op
p
µ(D)/c0.
Proof. Property (1) is immediate from the definition. For property (2), note that the unnormalized kernel ˜κ(x, y) :=
exp(−dϵ(Π(q(x)),Π(k(y)))/τ)satisfies˜κ≤1. Thus
R
D κ(x, y)dµ(x) =
R
D ˜κ(x, y)/Z(x)dµ(x)≤µ(D)/c 0.
28

Hyperbolic Neural Operator
Proposition F.3(A sufficient lower bound on the normalizer Z(x)).Assume there exist R <∞ and τmin >0 such that for
allx, y∈D,Π(q(x;a)),Π(k(y;a))∈B H(o, R)andτ≥τ min. Then
Z(x)≥µ(D) exp

− max{2R,arcosh(1 +ϵ)}
τmin

=: c0 for allx∈D.
In particular, it suffices to bound Euclidean feature norms: if ∥q(x;a)∥ ≤B q and ∥k(y;a)∥ ≤B k for all x, y, then one can
takeR= max{asinh(B q),asinh(B k)}.
Proof. If Π(q(x;a)),Π(k(y;a))∈B H(o, R), then the hyperbolic triangle inequality gives dH(Π(q(x;a)),Π(k(y;a)))≤
2R. By Proposition A.1 and monotonicity ofarcosh,
dϵ

Π(q(x;a)),Π(k(y;a))


= arcosh

max{1 +ϵ,−⟨Π(q(x;a)),Π(k(y;a))⟩ M }


≤max{2R,arcosh(1 +ϵ)}.
Hence exp(−dϵ/τ)≥exp(−max{2R,arcosh(1 +ϵ)}/τ min) for all y, and integrating over D yields the claimed bound
on Z(x). Finally, r(z) =d H(o,Π(z)) = asinh(∥z∥) implies Π(z)∈B H(o,asinh(∥z∥)) , so a uniform Euclidean norm
bound implies a uniform radius bound.
RemarkF.4 (How to enforce the assumptions in practice).The radius bound can be guaranteed by explicitly bounding ∥q∥
and ∥k∥ via bounded activations such as tanh, explicit norm clipping, or layer normalization followed by a bounded scale.
Ifτis learnable, a convenient way to ensureτ≥τ min is to parameterizeτ=τ min + softplus(ˆτ).
F.3. Bi-Lipschitz Residuals and Stability
For the main text, we analyze afixed-kernelsetting for the operator family {Ka} introduced in Section 4.2. We treat
the query and key maps as fixed for a given input field a. Under this conditioning, the kernel κa(x, y) is fixed and the
map v7→ K av is linear. In particular, it is ∥Ka∥-Lipschitz in any norm, and in L2(D) we have Lip(Ka) =∥K a∥L2→L2.
Combining this with the boundedness results from Appendix F.2 yields a sufficient uniform Lipschitz bound to invoke
Proposition F.6; the resulting condition is sufficient and can be conservative.
RemarkF.5 (Implementation (stabilization and temperature lower bound)).In our implementations, the per-head temperature
is clamped to τ∈[0.1,3.0] , and the arcosh argument is clamped to ≥1 +ϵ with ϵ∈ {10 −6,10 −4} to avoid the singular
derivative at 1. We also apply LayerNorm to the hidden state before Q/K projections. Tab. 7 confirms that the arcosh clamp
is never activated on benchmark evaluations of our best checkpoints.
Proposition F.6(Bi-Lipschitz stability of residual maps).Let (X,∥ · ∥) be a normed vector space and let F:X→X be
CF-Lipschitz:∥F(u)− F(v)∥ ≤C F ∥u−v∥. If|α|C F <1andT:=I+αF, then for allu, v∈X,
(1− |α|C F)∥u−v∥ ≤ ∥T(u)−T(v)∥ ≤(1 +|α|C F)∥u−v∥.
In particular,Tis injective and has a Lipschitz inverse on its range.
Proof.By the triangle inequality:
∥T(u)−T(v)∥=∥(u−v) +α(F(u)− F(v))∥
≤ ∥u−v∥+|α|∥F(u)− F(v)∥ ≤(1 +|α|C F)∥u−v∥,
∥T(u)−T(v)∥ ≥ ∥u−v∥ − |α|∥F(u)− F(v)∥ ≥(1− |α|C F)∥u−v∥.
This gives the claimed bi-Lipschitz bounds.
Application to Fixed-Kernel Layers.For each a, apply Proposition F.6 with X=L 2(D) and F=K a, so that
CF =∥K a∥L2→L2. Since Ka is linear, this condition is sharp up to constants and matches standard sufficient conditions for
invertible residual blocks (Behrmann et al., 2019).
G. Discretization Consistency
We analyze discretization consistency in thefixed-kernelsetting: for each input a, the feature maps q(·;a) and k(·;a) (hence
κ(·,·;a) ) are treated as fixed, and the discretization error is measured relative to the induced continuum operator. We first
prove that the discrete scalar-kernel operator Kκ,N converges to the continuum kernel operator Kκ as the mesh refines. The
full layer discretization is then KN :=W O ◦ Kκ,N ◦W V , which converges to K=W O ◦ Kκ ◦W V . This is crucial for
neural operators: it ensures the learned model does not depend on a particular discretization.
29

Hyperbolic Neural Operator
G.1. Setup and Notation
Let D⊂R d be a bounded domain with Lebesgue measure µ. Consider a sequence of discretizations {xj}N
j=1 ⊂D
with associated quadrature weights {wj}N
j=1. We assume the quadrature rule satisfies the Lipschitz accuracy condition
(Assumption G.3), i.e., for any Lipschitz integrandg:

Z
D
g(y)dµ(y)−
NX
j=1
wjg(xj)

≤C q hLip(g),(28)
where h:= max j supx∈Vj ∥x−x j∥ is the mesh size (maximum cell diameter), and Vj is the V oronoi cell ofxj. We equip
vectorsu∈(R dv)N with the weighted norm∥u∥ 2
ℓ2w
:=PN
j=1 wj∥uj∥2
2.
The continuum operator is:
(Kκv)(x) =
Z
D
κ(x, y)v(y)dµ(y),
and we define the sampling operator PN by (PN v)j :=v(x j). The Nyström (quadrature) discretization of the continuum
kernel is the mapK cont
κ,N : (Rdv)N →(R dv)N given by
(Kcont
κ,N u)i :=
NX
j=1
wj κ(xi, xj)u j, u∈(R dv)N .
The implemented discrete scalar operator Kκ,N uses the discretely normalized kernel κN from Eq. 29 and applies it with
quadrature weights,
(Kκ,N u)i :=
NX
j=1
wj κN(xi, xj)u j =
NX
j=1
αij uj, α ij =w j κN(xi, xj).
The full layer discretization isK N :=W O ◦ Kκ,N ◦W V .
We prove discretization consistency.
G.2. Kernel Regularity
We first establish regularity of the hyperbolic attention kernel.
Lemma G.1(Lipschitzness of the fixed-kernel hyperbolic attention kernel).Let κ(x, y) = ˜κ(x, y)/Z(x)with ˜κ(x, y) =
exp(−dϵ(Π(q(x)),Π(k(y)))/τ) and dϵ(p, q) = arcosh(max{−⟨p, q⟩M ,1 +ϵ}) . Assume that q and k are Lipschitz with
constants Lq and Lk, that Π(q(D)) and Π(k(D)) lie in a hyperbolic ball BH(o, R), and that Z(x)≥c 0 >0 and
τ≥τ min >0hold uniformly. LetL qk := max{Lq, Lk}. Then for each fixedx, the mapy7→˜κ(x, y)is globally Lipschitz
onDwith
Lipy(˜κ(x,·))≤C R · Lqk
τmin
· 1p
ϵ(2 +ϵ)
.
Consequently,y7→κ(x, y)is globally Lipschitz onDwith
Lipy(κ(x,·))≤C R · Lqk
τmin c0
· 1p
ϵ(2 +ϵ)
,
and forϵ >0the kernel is piecewiseC ∞ away from the clamping boundary{−⟨p, q⟩ M = 1 +ϵ}and globally Lipschitz.
Proof.Proof outline.Defineg(t) := arcosh(max{t,1 +ϵ})forϵ >0. Since d
dt arcosh(t) = 1√
t2−1 fort >1, one has
|g′(t)| ≤ 1p
(1 +ϵ) 2 −1
= 1p
ϵ(2 +ϵ)
=:L g,
and the clamping makesggloballyL g-Lipschitz onR.
30

Hyperbolic Neural Operator
Next, the Lorentz lift Π(z) = (
p
1 +∥z∥ 2, z) satisfies that the scalar map z7→
p
1 +∥z∥ 2 has gradient norm ≤1 , hence
it is1-Lipschitz. Therefore, for allz 1, z2 ∈R m,
∥Π(z1)−Π(z 2)∥2
2 =
p
1 +∥z 1∥2 −
p
1 +∥z 2∥2
2
+∥z 1 −z 2∥2
2 ≤2∥z 1 −z 2∥2
2,
soΠisL Π :=
√
2-Lipschitz.
(i) Stabilized distance is Lipschitz on bounded sets.Recall dϵ(p, q) =g(−⟨p, q⟩ M), and |⟨a, b⟩M | ≤ ∥a∥ 2∥b∥2. Fix
x∈D and write px := Π(q(x;a)) . If px and my := Π(k(y;a)) remain in a hyperbolic ball BH(o, R), then their Euclidean
norms are uniformly bounded; for instance
∥w∥2 ≤coshR+ sinhR≤e R =: BR for allw∈B H(o, R).
For anyy 1, y2 ∈D,
dϵ(px, my1)−d ϵ(px, my2)
 ≤L g
⟨px, my2 −m y1 ⟩M
 ≤L g ∥px∥2 ∥my2 −m y1 ∥2 ≤L g BR ∥my2 −m y1 ∥2.
Using the Lipschitzness ofk(·;a)andΠgives
∥my2 −m y1 ∥2 ≤L Π Lk ∥y2 −y 1∥2,
and similarly with Lq when varying x. Thus, for each fixed x, the map y7→d ϵ(Π(q(x;a)),Π(k(y;a))) is Lipschitz with
constant≤L g BR LΠ Lqk.
(ii) Kernel Lipschitzness and normalization.For τ≥τ min >0 , the map t7→e −t/τ has derivative bounded by 1/τmin on
[0,∞), and˜κ(x, y)∈(0,1]. Hence
Lipy(˜κ(x,·))≤ 1
τmin
Lipy

dϵ(Π(q(x;a)),Π(k(·;a)))


≤ Lg BR LΠ Lqk
τmin
.
Finally, κ(x, y) = ˜κ(x, y)/Z(x)with Z(x)≥c 0 >0 by assumption, and for fixed x the denominator is constant w.r.t. y.
Therefore
Lipy(κ(x,·))≤ 1
c0
Lipy(˜κ(x,·))≤ Lg BR LΠ Lqk
τmin c0
.
This yields the stated bound (absorbing BRLΠ into the radius-dependent constant CR) and also shows Lipschitzness for
ϵ >0under the bounded-radius assumption.
Lemma G.2(Feature Map Lipschitz Bound).Assume, for a fixed input field a and layer state, that the feature maps
satisfy q(·;a), k(·;a)∈W 1,∞(D) with constants Lq and Lk. This condition holds, for example, when the input fields
and positional features entering the Q/K maps are Lipschitz and the finite-dimensional projection networks have bounded
Lipschitz constants. With weight normalization or spectral normalization, the network contribution to these constants can be
controlled.
G.3. Quadrature Error Analysis
For normalized kernels, the discrete operator uses softmax normalization:
κN(xi, xj) := exp(−dϵ(Π(q(xi)),Π(k(x j)))/τ)PN
ℓ=1 wℓ exp(−dϵ(Π(q(xi)),Π(k(x ℓ)))/τ)
.(29)
This differs from the continuum kernel which uses integral normalization. In our experiments, we take uniform quadrature
weights wj =µ(D)/N (equivalently wj = 1 up to a global scaling), so the implemented coefficients αij =w j κN(xi, xj)
reduce to standard softmax attention weights. We keep the weighted form for general nonuniform discretizations.
Assumption G.3(Lipschitz quadrature accuracy).Assumew j ≥0and PN
j=1 wj =µ(D). There existsC q >0such that
for any Lipschitz functiongonD,

Z
D
g(y)dµ(y)−
NX
j=1
wjg(xj)

≤C q hLip(g).
31

Hyperbolic Neural Operator
RemarkG.4 (Relation to W1).Assumption G.3 holds, for instance, when the normalized measures ¯µ:=µ/µ(D) and
¯µN := 1
µ(D)
PN
j=1 wjδxj satisfy W1(¯µ,¯µN)≤Ch for some constant C independent of N, by the Kantorovich–Rubinstein
duality (Villani, 2008). This includes quasi-uniform meshes, and can be interpreted as a mesh-quality condition in point-cloud
settings with uniform weights.
Theorem G.5(Discretization consistency under Lipschitz kernels).Assume Lemma G.1 and Assumption G.3. If v is
bounded and Lipschitz onD, then
∥Kκ,N PN v−P N Kκv∥ℓ2w ≤C h∥v∥ W 1,∞(D),
where C depends on Cq, the uniform bound on κ, and supx Lipy(κ(x,·)) . All norms are understood in the vector-valued
sense using the Euclidean norm in the feature space.
Consequently, the full layer discretization satisfies
∥KN PN v−P N Kv∥ℓ2w ≤ ∥W O∥op∥WV ∥op C h∥v∥ W 1,∞(D).
Proof.We decompose the error at sample points into two terms:
Kκ,N PN v−P N Kκv=

Kκ,N PN v− K cont
κ,N PN v


| {z }
(I) Normalization gap
+

Kcont
κ,N PN v−P N Kκv


| {z }
(II) Quadrature
.
(I) Normalization gap.Let ˜κi(y) := exp(−d ϵ(Π(q(xi)),Π(k(y)))/τ) , Z(x i) :=
R
D ˜κi(y)dµ(y) , and ZN(xi) :=PN
j=1 wj˜κi(xj). Thenκ(x i, xj) = ˜κi(xj)/Z(xi)andκ N(xi, xj) = ˜κi(xj)/ZN(xi), hence
(Kκ,N PN v)i −(K cont
κ,N PN v)i
 ≤ ∥v∥ ∞

1
ZN(xi) − 1
Z(x i)

NX
j=1
wj˜κi(xj)
=∥v∥ ∞
|Z(x i)−Z N(xi)|
Z(x i) ≤ ∥v∥∞
c0
|Z(x i)−Z N(xi)|.
Applying Assumption G.3 to y7→˜κ i(y) gives |Z(x i)−Z N(xi)| ≤C qhLip(˜κi). By Lemma G.1, Lip(˜κi) is uniformly
bounded (and does not incur an additional 1/c0 factor). SincePN
j=1 wj =µ(D)>0 and ˜κi(xj)>0 , we have ZN(xi)>0
andκ N is well-defined. Thus (I) is bounded byCh∥v∥ ∞ for a constantCindependent ofN.
(II) Quadrature.Fixiand defineg i(y) :=κ(x i, y)v(y). By Assumption G.3,
(Kcont
κ,N PN v)i −(P N Kκv)i
 =

NX
j=1
wjgi(xj)−
Z
D
gi(y)dµ(y)

≤C qhLip(g i).
Moreover, Lip(gi)≤ ∥κ(x i,·)∥ ∞ Lip(v) +∥v∥ ∞ Lipy(κ(xi,·)) . Using Lemma G.1 to bound supi Lipy(κ(xi,·)) and the
uniform bound ∥κ∥∞ ≤1/c 0, summing over quadrature points with weights wi and usingP
i wi =µ(D) , yields a bound
≤C
p
µ(D)h∥v∥ W 1,∞ inℓ 2
w for a constantCindependent ofN.
Combining (I) and (II) and using∥v∥ ∞ ≤ ∥v∥ W 1,∞(D) yields the stated estimate.
RemarkG.6 (Discrete normalization stability).Since PN
j=1 wj =µ(D)>0 and ˜κi(xj)>0 , one always has ZN(xi)>0 .
Moreover, if h is small enough such that Cqhsup i Lip(˜κi)≤c 0/2, then ZN(xi)≥c 0/2 for all i, which improves
numerical stability of the discrete normalization.
Corollary G.7(Strong convergence / discretization invariance).Under the assumptions of Theorem G.5, for any sequence
of discretizations with mesh sizeh→0one has, for eachv∈W 1,∞(D;R dv),
∥Kκ,N PN v−P N Kκv∥ℓ2w →0,∥K N PN v−P N Kv∥ℓ2w →0.
Equivalently, the difference operators converge strongly and haveO(h) error as maps W 1,∞(D;R dv)→ℓ 2
w, matching the
continuum-referenced discretization consistency commonly required in neural operators (Kovachki et al., 2023).
Proof.Immediate from Theorem G.5.
32

Hyperbolic Neural Operator
G.3.1. PRACTICALIMPLICATIONS
RemarkG.8 (Convergence Rate).For typical meshes withN∼h −d points, Theorem G.5 implies:
∥Kκ,N PN v−P N Kκv∥ℓ2w =O(N −1/d).
RemarkG.9 (Stability of Hyperbolic Distance).The arcosh function has a singularity at z= 1 (when dH = 0). In our
implementation, we clamp the argument to z≥1 +ϵ with a small ϵ (typically 10−6 or 10−4 depending on the variant),
which:
1. Ensures the derivative1/
√
z2 −1is bounded byO(1/ √ϵ).
2. Introduces a bias of orderO( √ϵ)in the distance, which is negligible.
This numerical stabilization does not affect the convergence rate.
RemarkG.10 (Temperature Dependence).The Lipschitz constant in Lemma G.1 depends on 1/τ. Lower temperatures
τ→0 sharpen the attention but increase the kernel’s Lipschitz constant, potentially degrading the discretization error
constant. In our implementations, we clamp the learned temperatures toτ∈[0.1,3.0]for numerical stability.
H. Toy Multiscale Tree-Kernel Fitting: Quantitative Scaling
We report quantitative errors for the toy tree-kernel fitting experiment summarized in the main text.
Setup.This toy is a controlledgraph-domainoperator-learning problem (binary tree), not a physical PDE benchmark. We
use complete binary trees of depthL∈ {5,6,7,8,9,10} with N=2 L leaves. Let dtree(i, j) denote the shortest-path distance
between leaves i and j. Equivalently, if ℓ(i, j)∈ {0, . . . , L} is the depth of the lowest common ancestor (LCA) of leaves
i, j, then dtree(i, j) = 2(L−ℓ(i, j)) , which induces L+1 interaction scales by LCA depth. We define the ground-truth
row-stochastic kernel
Ktrue(i, j)∝exp

−γ eff dtree(i, j)


, γ eff =γ/L,
with γ= 1 . We fit Kpred from learned leaf embeddings using either dot-product or hyperbolic-distance attention (both with
embedding dimension d=4 and learnable temperature τ). For hyperbolic attention, we use a Poincaré-ball distance with
learnable curvature c (bounded by a sigmoid parameterization for stability). We optimize a mixed objective combining
row-wise kernel divergenceKL(Ktrue∥Kpred) and an operator regression loss on random sourcess∼ N(0, I) , using AdamW
(lr10 −2, weight decay10 −4) for 2000 steps with batch size 64, and keep the optimization budget fixed acrossN.
Metrics.We report the relative kernel error ∥Kpred −K true∥F /∥Ktrue∥F , the relative operator error Es

∥Kpreds−
Ktrues∥2/∥Ktrues∥2

withs∼ N(0, I), and the row-wise KL divergenceKL(K true∥Kpred).
32 64 128 256 512 1024
N
0.1
0.2
0.3Error
Rel. Frobenius
32 64 128 256 512 1024
N
0.1
0.2
0.3
Rel. 2
32 64 128 256 512 1024
N
0.00
0.01
0.02
0.03
Row KL
Hyperbolic Dot-product
Figure 10.Toy tree-kernel fitting: quantitative scaling with leaf countN.
33

Hyperbolic Neural Operator
Table 8.Numerical values for Fig. 10.
N REL. FROBENIUS REL.ℓ 2 ROWKL
HYPDOT HYPDOT HYPDOT
32 0.040 0.298 0.039 0.273 0.0006 0.0361
64 0.050 0.321 0.050 0.322 0.0008 0.0363
128 0.057 0.317 0.053 0.312 0.0010 0.0323
256 0.087 0.298 0.082 0.288 0.0026 0.0269
512 0.108 0.271 0.110 0.255 0.0032 0.0219
1024 0.096 0.243 0.106 0.261 0.0023 0.0176
PDE Viewpoint and Properties.Although this experiment is phrased as kernel fitting, the mapping u=K trees can be
viewed as a discrete nonlocal diffusion step on a hierarchical domain (the leaves of a binary tree). Since Ktree is nonnegative
and row-stochastic, each output entry is a convex combination of inputs, so it obeys a discrete maximum principle:
min
j
sj ≤(K trees)i ≤max
j
sj,∥K trees∥∞ ≤ ∥s∥ ∞, K tree1=1.
This makes the operator stable and mass-preserving on constant modes, consistent with diffusion-type PDE solution operators.
Our exponential distance kernel Kij ∝exp(−γ effdtree(i, j)) mirrors the off-diagonal decay typical of Green’s/heat kernels
of elliptic operators on tree-like geometries. Scaling γeff =γ/L keeps the effective interaction range comparable as depth
grows (sinced tree scales linearly inL), avoiding trivial near-diagonal or near-uniform limits at largeN.
Connection to a Graph Elliptic PDE.Distance-decaying kernels on trees also arise from diffusion/elliptic operators on
graphs. For example, consider the (tree-)graph elliptic problem
(I+αL+βL 2)u=s,(30)
where L=D−A is the graph Laplacian of the full binary tree. For α, β≥0 the operator is symmetric positive definite, so
the solution exists and is unique, and it acts as a spectral low-pass filter: if Lϕk =λ kϕk, then ⟨u, ϕk⟩=g(λ k)⟨s, ϕk⟩ with
g(λ) = 1/(1 +αλ+βλ 2). The resulting smoothing couples leaves across multiple scales induced by shared ancestors,
which is precisely the hierarchical structure the toy is designed to expose.
I. Training and Implementation Details
Model Variants.HNO uses the same hyperbolic-distance Gibbs kernel across all datasets, with tokenization chosen to
match the underlying discretization. For point clouds (Elasticity), we use a learned latent set of M tokens and perform
latent↔point cross-attention to encode/decode, with hyperbolic self-attention operating on the latent tokens. For regular
grids and structured meshes (Navier–Stokes, Darcy, Airfoil, Pipe), we tokenize the field into non-overlapping 2D patches
via a strided convolutional embedding; attention is computed between patch tokens (with a 2D relative position bias), and
predictions are mapped back to the original resolution with an unpatching operator. For the spatiotemporal Plasticity setting,
we keep the same patch tokenization and hyperbolic attention, but use a stronger local feature extractor to better preserve
fine-scale dynamics.
Computational Complexity.Let N be the number of input discretization nodes and let M be the number of tokens
processed by the interaction core. Dense attention over all nodes incurs O(N 2) pairwise interactions per layer. For point
clouds, we use a latent set ofM tokens with cross-attention for encoding and decoding, which costsO(N M), and hyperbolic
self-attention on the latent set, which costs O(M 2). This yields an overall interaction cost O(N M+M 2) per layer. For
regular grids and structured meshes, patchification produces M patch tokens and self-attention runs in the patch space with
cost O(M 2), where M≪N at fixed patch size. Replacing dot-product logits by stabilized hyperbolic-distance logits does
not change these asymptotic costs, since computing Minkowski inner products andarcoshis constant-time per token pair.
Training Protocol.Unless otherwise noted, we train PDEBench benchmarks for 500 epochs; Elasticity uses 300 epochs
and Airfoil uses 800 epochs. For PDEBench, we follow the standard split of 1000 training and 200 test samples and report
mean relative ℓ2 error (averaged over three runs). We use AdamW with weight decay 10−5 and apply gradient clipping
34

Hyperbolic Neural Operator
Table 9.Best hyperparameters for HNO on PDEBench. Here d is hidden width, L is depth, H is the number of heads, and dH is the
hyperbolic Q/K dimension. Tokenization is either a learned latent set (Mtokens) or patchification with patch sizeps.
DATASET TOKENS d L H dH LR EPOCHS BATCH CLIP SCHED.
ELASTICITY M=96 384 6 8 16 5×10 −4 300 4 – COSINE
NAVIER–STOKES ps=4 132 4 4 16 1×10 −3 500 8 0.5 ONECYCLE
DARCY ps=5 96 4 4 16 1×10 −3 500 4 0.1 ONECYCLE
PLASTICITY ps=3 128 5 8 16 8×10 −4 500 16 0.1 ONECYCLE
AIRFOIL ps=5 56 5 4 16 1×10 −3 800 4 0.5 ONECYCLE
PIPE ps=3 132 4 4 16 1×10 −3 500 4 1.0 ONECYCLE
Table 10.Top-3hyperparameter configurations per PDEBench dataset (lower Rel.ℓ 2 is better; best per dataset in bold).
DATASET TOKENS d L H dH LR PARAMS REL.ℓ 2
ELASTICITY M=96 384 6 8 16 5×10 −4 15.37M 0.00368971
ELASTICITY M=64 320 8 8 16 1×10 −3 14.23M 0.00374535
ELASTICITY M=64 256 10 8 16 1×10 −3 11.47M 0.00410951
NAVIER–STOKES PS=4 132 4 4 16 1×10 −3 1.19M 0.0675513
NAVIER–STOKES PS=4 128 4 4 16 1.2×10 −3 1.12M 0.0679059
NAVIER–STOKES PS=4 128 4 4 16 1×10 −3 1.12M 0.0689504
DARCY PS=5 96 4 4 16 1×10 −3 0.82M 0.00446
DARCY PS=5 96 5 4 16 1×10 −3 0.91M 0.00446
DARCY PS=5 104 4 4 16 1×10 −3 0.91M 0.004462
PLASTICITY PS=3 128 5 8 16 8×10 −4 2.27M 0.000903203
PLASTICITY PS=5 128 4 8 16 8×10 −4 2.38M 0.00107384
PLASTICITY PS=3 128 4 8 16 5×10 −4 1.89M 0.0011549
AIRFOIL PS=5 56 5 4 16 1×10 −3 0.36M 0.00508543
AIRFOIL PS=5 96 6 4 16 5×10 −4 1.17M 0.00510241
AIRFOIL PS=5 112 5 4 16 5×10 −4 1.23M 0.00515098
PIPE PS=3 132 4 4 16 1×10 −3 5.10M 0.00273442
PIPE PS=3 128 3 4 16 1×10 −3 3.86M 0.00284324
PIPE PS=3 132 4 4 16 1.2×10 −3 5.10M 0.0029343
Discussion (Elasticity).In our Elasticity sweeps, the best-performing configurations use larger hidden widths (and thus larger parameter counts), which we attribute to the harder
encoding problem induced by irregular, nonuniform point sampling. This is orthogonal to overly smooth tokenizations that can suppress high-frequency modes; we believe
designing geometry-aware tokenizers/encoders that better preserve fine-scale content on irregular 2D/3D point sets is a valuable direction, but we do not explore it in this work.
when needed for stability. Elasticity uses cosine learning-rate annealing, while the other PDEBench datasets use a OneCycle
learning-rate schedule (stepped per batch).
Preprocessing and Normalization.For Darcy/Pipe/Elasticity, we standardize targets (and inputs where applicable) using
dataset statistics and decode predictions back to the original scale before reporting relative ℓ2. For Airfoil, we do not
normalize the output and (for our best model) keep the input coordinates unnormalized due to the native physical-coordinate
range.
PDEBench Benchmarks.Tab. 9 summarizes the best hyperparameters used for HNO on the six PDEBench benchmarks.
Hyperparameter Search Robustness (Top-3).To demonstrate robustness and reduce concerns about configuration
cherry-picking, we report the top-3 configurations (sorted by test relativeℓ2) from our lightweight sweep on each benchmark,
using the same data split and evaluation metric as in the main results. Each row lists the tokenization parameter (M or patch
size ps), architecture hyperparameters (d, L, H, dH), learning rate, parameter count, and the corresponding test error. These
rows summarize lightweight sweep configurations used to select architectures; final main-table results may use the selected
configuration under the main evaluation protocol.
Large-Scale Benchmarks.We implement HNO on large-scale unstructured meshes by replacing dot-product logits in
Transolver++ with stabilized hyperbolic-distance logits, keeping the same backbone width and depth. We use d=256,
L=8, H=8, 64 slices, and dH=16, trained with Adam and a OneCycle learning-rate schedule and batch size 1. Tab. 11
summarizes the dataset-specific settings.
35

Hyperbolic Neural Operator
Table 11.Large-scale unstructured-mesh benchmarks: training setup.
DATASET #NODES LR EPOCHS GRAPH SUBSAMPLE NOTES
SHAPENETCAR 32,186 1×10 −3 500 RADIUS(r=0.2) – LOSS IS VELOCITYMSEPLUS WEIGHTED SURFACE-PRESSUREMSE
(λ=0.5).
AIRFRANS 32,000 1×10 −3 500 RADIUS(r=0.05) 32K PER-EPOCH RESAMPLING WITH MAX NEIGHBORS64;WEIGHTEDMSE
(SURFACE+VOLUME).
Table 13.Mixed-precision stability check on Darcy (85×85, batch 4).
PRECISION MODE STEPS NANS
FP16 TRAIN(AMP + GRADSCALER) 200 NO
BF16 FORWARD-ONLY 200 NO
Visualization of Large-Scale Error Maps.For AirfRANS (2D), we load the ground-truth internal mesh and model
predictions (VTK .vtu), compute signed errors in velocity magnitude and pressure, and render them as filled-contour maps
using the original mesh connectivity (triangulating quads when needed). Color limits are set per sample using symmetric
percentile clipping and shared across models for fair visual comparison; the airfoil boundary is overlaid for context. For
ShapeNet Car (3D), we render a semi-transparent outer-domain surface colored by signed speed residual (∥Upred∥ − ∥Ugt∥)
together with a surface pressure error map (ppred −p gt); meshes are aligned to a consistent upright orientation across cases.
J. Additional Experiments
J.1. Ablations and Sensitivity
Darcy sensitivity.We sweep hidden dimension, the number of attention layers, and hyperbolic Q/K dimension while
keeping other settings fixed (Fig. 8). Tab. 3 reports the official baseline configurations, which are not parameter-matched to
HNO; for a parameter-matched efficiency microbenchmark, see Tab. 12.
Table 12.Parameter-matched efficiency microbenchmark on Darcy (A6000 48GB, 85×85 grid, batch 4). “Train” is extrapolated from the
measured per-step time to 500 epochs with 1000 training samples.
METHOD PARAMS
VRAM
(GB)↓
TRAIN
(H)↓
INFER
(MS/B)↓
HNO (OURS) 0.82M 0.22 1.18 5.69
TRANSOLVER++ (OFFICIAL) 2.84M 3.07 4.03 54.79
TRANSOLVER++ (PARAM-MATCHED) 0.82M 1.53 1.97 50.84
Parameter-matched efficiency on Darcy.To isolate architectural efficiency from model size, we benchmark a parameter-
matched Transolver++ configuration on Darcy by reducing hidden width and depth while keeping the slicing mechanism
unchanged. We report peak allocated GPU memory during a single training step, along with per-batch inference latency. All
measurements use an NVIDIA A6000 (48GB), a85×85grid (downsample 5), and batch size 4.
J.2. Numerical Stability and Runtime Breakdown
Mixed-precision stability (Darcy).We run a short stability check on Darcy ( 85×85, batch 4). FP16 uses AMP with
GradScaler for 200 optimization steps and produces no NaNs. Our current PyTorch build does not support BF16 backprop
throughrsqrt, so we report a BF16 forward-only check.
Runtime breakdown (Darcy).We profile a forward pass on Darcy (85×85, batch 4) and report a coarse operator-level time
breakdown. Hyperbolic distance nonlinearity (arcosh) accounts for<1%of self CUDA time.
J.3. Spectral Error Analysis on Darcy
We analyze the frequency content of prediction errors on Darcy (85×85) by computing a radial average of the relative error
spectrum over 200 test samples: for each radial bin b, we report
qP
∥k∥∈b|ˆe(k)|2 /P
∥k∥∈b|ˆu(k)|2, where e=u pred −u gt.
36

Hyperbolic Neural Operator
Table 14.Runtime breakdown on Darcy ( 85×85, batch 4). Percentages are relative to total self CUDA time from a CUDA profiler trace;
latency is measured by CUDA events.
METHOD
LATENCY
(MS/B)↓
arcosh
(%)
MATMUL+EINSUM
(%)
CONV
(%)
HNO (HYPERBOLIC) 5.87 0.8 15.7 25.0
−HYPERBOLIC(EUCLID) 5.18 0.0 16.7 22.1
0 10 20 30 40 50
Radial frequency |k|
10 2
10 1
100
Relative error spectrum
FNO
FixedTau
Euclid
HNO
Figure 11.Darcy error spectrum (log scale). HNO yields uniformly lower spectral error than Euclidean-distance and fixed-temperature
ablations, and does not exhibit increased high-frequency error relative to an FNO baseline.
J.4. Causal Radius Interventions
Motivation.Fig. 9 shows that learned hyperbolic radii correlate with attention locality. To test whether radius is acausal
control signal rather than a byproduct of training, we perform radius interventions that preserve the Q/K directions but
disrupt their norms.
Interventions.Let qi ∈R dH and kj ∈R dH denote Euclidean features that are lifted to the Lorentz model. Write qi =∥q i∥¯qi
and kj =∥k j∥¯kj with ∥¯qi∥=∥ ¯kj∥= 1 , so the radius is r(qi) = asinh(∥qi∥). We consider two interventions applied at
inference time: (i)Fixed norm:set all ∥qi∥ (and all ∥kj∥) to their within-head mean while keeping ¯qi, ¯kj unchanged; (ii)
Shuffled norm:randomly permute the set of {∥qi∥}i (and {∥kj∥}j) within each head, again keeping directions unchanged.
Both interventions keep angles in feature space and the value pipeline fixed, while removing meaningful radius stratification.
Results (Darcy).Tab. 15 reports the mean relative ℓ2 error on Darcy (85×85). Both interventions substantially degrade
accuracy, indicating that learned radii are important for this checkpoint’s near–far routing behavior.
Table 15.Causal radius interventions on Darcy (85×85).
METHOD REL.ℓ 2 ↓
HNO (BASELINE) 0.00446
FIXED NORM(Q/K) 0.03269
SHUFFLED NORM(Q/K) 0.03328
J.5. Learned Kernel Decay
Assumption D.2 posits that the far-field kernel is smooth and decays away from the near field. As an empirical sanity check,
we visualize distance-dependent decay in learned attention weights on Elasticity. For a representative global token (small
hyperbolic radius) and a representative local token (large hyperbolic radius), we compute each token’s physical barycenterb
and plot log10 attention weights against the physical distance ∥x−b∥ 2. Local tokens exhibit a much faster decay, while
global tokens remain broadly connected.
37

Hyperbolic Neural Operator
0.2 0.3 0.4 0.5 0.6 0.7 0.8
Distance to token barycenter x b 2
4.0
3.5
3.0
2.5
2.0
1.5
log10  attention weight
Global token (latent 14, d=0.19)
Binned median
0.0 0.2 0.4 0.6 0.8 1.0
Distance to token barycenter x b 2
Local token (latent 55, d=0.90)
Binned median
Figure 12.Elasticity: distance-dependent decay of learned attention weights for representative global and local tokens. The curve shows
the binned median oflog 10 attention weights over distance to the token-specific barycenter.
J.6. Discretization Robustness
Quadrature weights in the discrete operator.Our Nyström discretization uses nodes {xj}N
j=1 with nonnegative weights
{wj}N
j=1 (Eq. (11)). On PDEBench, we use uniform weights wj =µ(D)/N (equivalently wj = 1 up to a global scaling),
so the implemented coefficients reduce to standard softmax attention weights.
Scope of the consistency result.The discretization consistency theorem (Appendix G) assumes a Lipschitz-accurate
quadrature rule (Assumption G.3), which can be interpreted as a mesh-quality condition. This covers quasi-uniform grids
and well-behaved point clouds whose empirical measures converge to the underlying continuum measure in W1 at rate
O(h). The theorem does not apply to arbitrary resampling or severe subsampling that violates this quadrature accuracy
condition. Empirically, our reported benchmarks evaluate each dataset on its native discretization. To provide a concrete
discretization-density sanity check, Appendix J.6.1 additionally reports an AirfRANS node subsampling scan that varies
the available mesh nodes at test time. We leave robustness to aggressive remeshing schemes or distribution-shifted point
sampling as future work.
J.6.1. AIRFRANS DISCRETIZATIONDENSITYSCAN
Setup.We evaluate Transolver-series models on AirfRANS by varying the number of nodes available at test time. For
a representative test case, we uniformly subsample the mesh nodes and report results at representative node budgets
N∈ {2k,32k,128k} . We report mean relative ℓ2 error in physical units, peak GPU memory, and end-to-end forward
latency on an RTX A6000 48GB.
Table 16.AirfRANS subsampling scan on a representative test case. Each entry reports mean relative ℓ2 / peak VRAM (GB) / forward
latency (ms) at the given node budgetN.
N TRANSOLVERTRANSOLVER++ HNO
2K 0.509 / 0.28 / 126.0 0.334 / 0.29 / 44.9 0.165 / 0.31 / 76.4
32K 0.478 / 3.49 / 148.9 0.237 / 3.75 / 171.1 0.106 / 3.98 / 186.0
128K 0.469 / 13.65 / 533.6 0.235 / 14.72 / 535.2 0.105 / 15.59 / 515.7
K. Additional Visualizations
K.1. Qualitative Results on PDEBench
We provide additional qualitative visualizations for all six PDEBench datasets using the HNO best checkpoints. Each row
corresponds to one test case; columns show input (when available), ground truth (GT), HNO prediction, and error (prediction
minus GT).
38

Hyperbolic Neural Operator
Input
 Ground Truth
 H/glyph1197O Pred.
 Error
Figure 13.Darcy qualitative results (page 1/2). Columns show Input / GT / HNO prediction / Error (prediction minus GT).
39

Hyperbolic Neural Operator
Input
 Ground Truth
 H/glyph1197O Pred.
 Error
Figure 14.Darcy qualitative results (page 2/2). Columns show Input / GT / HNO prediction / Error (prediction minus GT).
40

Hyperbolic Neural Operator
Input
 Ground Truth
 H/glyph1197O Pred.
 Error
Figure 15.Airfoil qualitative results (page 1/2). Columns show Input / Ground Truth / HNO prediction / Error (prediction minus GT).
41

Hyperbolic Neural Operator
Input
 Ground Truth
 H/glyph1197O Pred.
 Error
Figure 16.Airfoil qualitative results (page 2/2). Columns show Input / Ground Truth / HNO prediction / Error (prediction minus GT).
42

Hyperbolic Neural Operator
Input
 Ground Truth
 H/glyph1197O Pred.
 Error
Figure 17.Pipe qualitative results (page 1/2). Columns show Input / Ground Truth / HNO prediction / Error (prediction minus GT).
43

Hyperbolic Neural Operator
Input
 Ground Truth
 H/glyph1197O Pred.
 Error
Figure 18.Pipe qualitative results (page 2/2). Columns show Input / Ground Truth / HNO prediction / Error (prediction minus GT).
44

Hyperbolic Neural Operator
Ground Truth
 H/glyph1197O Pred.
 Error
Figure 19.Elasticity qualitative results (page 1/2). Columns show Ground Truth / HNO prediction / Error (prediction minus GT).
45

Hyperbolic Neural Operator
Ground Truth
 H/glyph1197O Pred.
 Error
Figure 20.Elasticity qualitative results (page 2/2). Columns show Ground Truth / HNO prediction / Error (prediction minus GT).
46

Hyperbolic Neural Operator
Ground Truth
 H/glyph1197O Pred.
 Error
Figure 21.Plasticity qualitative results (page 1/1). Columns show Ground Truth / HNO prediction / Error (prediction minus GT).
47

Hyperbolic Neural Operator
Ground Truth
 H/glyph1197O Pred.
 Error
Figure 22.Navier–Stokes qualitative results (page 1/2). Columns show Ground Truth / HNO prediction / Error (prediction minus GT).
48

Hyperbolic Neural Operator
Ground Truth
 H/glyph1197O Pred.
 Error
Figure 23.Navier–Stokes qualitative results (page 2/2). Columns show Ground Truth / HNO prediction / Error (prediction minus GT).
49

Hyperbolic Neural Operator
K.2. Direct Attention Map Visualization
Figures 24–31 visualize representative global (small radius) and local (large radius) query tokens on Elasticity using cached
attention weights. These examples support the qualitative trend that larger-radius tokens attend more sharply: their heatmaps
concentrate around a tight neighborhood and their distance-decay curves drop faster. At the same time, the separation is not
a hard global–local dichotomy. First, locality is acontinuousproperty: tokens with intermediate radii can exhibit partially
overlapping receptive fields, and even “global” tokens can be anisotropic (e.g., emphasizing particular boundary/geometry
regions). Second, overlap is expected when the selected queries are spatially close or when the solution requires both
near-field fidelity and weak far-field conditioning (e.g., to satisfy global constraints), which leaves nonzero long-range tails
after normalization. Overall, the hyperbolic kernel induces alearnable routing biasrather than an explicit partition: radius
interacts with angular alignment and temperature, so the effective interaction range adapts smoothly to the input and the
local geometry.
50

Hyperbolic Neural Operator
10 4
10 3
10 2
Attention weight (log scale)
Global token (latent 14, d=0.19, RF=0.38)
 Local token (latent 55, d=0.90, RF=0.20)
0.0 0.2 0.4 0.6 0.8 1.0
10 4
10 3
10 2
 Global (d=0.19)
Local (d=0.90)
Global token (latent 5, d=0.29, RF=0.40)
 Local token (latent 22, d=0.88, RF=0.22)
0.0 0.2 0.4 0.6 0.8
10 4
10 3
10 2
 Global (d=0.29)
Local (d=0.88)
Global token (latent 56, d=0.45, RF=0.41)
 Local token (latent 50, d=0.87, RF=0.22)
0.0 0.2 0.4 0.6 0.8
10 3
10 2
 Global (d=0.45)
Local (d=0.87)
Global token (latent 20, d=0.49, RF=0.39)
 Local token (latent 27, d=0.87, RF=0.21)
0.0 0.2 0.4 0.6 0.8
Physical distance to barycenter x b 2
10 3
10 2
 Global (d=0.49)
Local (d=0.87)
Mean attention weight
Elasticity attention maps: additional global vs local examples
Figure 24.Additional attention-map examples on Elasticity, page 1. Each row shows a global token (left), a local token (middle), and their
binned mean attention decay versus distance (right).
51

Hyperbolic Neural Operator
10 4
10 3
10 2
Attention weight (log scale)
Global token (latent 41, d=0.55, RF=0.34)
 Local token (latent 57, d=0.87, RF=0.18)
0.0 0.2 0.4 0.6 0.8 1.0
10 4
10 3
10 2
Global (d=0.55)
Local (d=0.87)
Global token (latent 40, d=0.56, RF=0.31)
 Local token (latent 18, d=0.87, RF=0.25)
0.2 0.4 0.6 0.8
10 3
10 2
 Global (d=0.56)
Local (d=0.87)
Global token (latent 30, d=0.57, RF=0.38)
 Local token (latent 31, d=0.87, RF=0.20)
0.0 0.2 0.4 0.6 0.8 1.0
10 4
10 3
10 2
 Global (d=0.57)
Local (d=0.87)
Global token (latent 4, d=0.59, RF=0.26)
 Local token (latent 11, d=0.87, RF=0.23)
0.0 0.2 0.4 0.6 0.8 1.0
Physical distance to barycenter x b 2
10 4
10 3
Global (d=0.59)
Local (d=0.87)
Mean attention weight
Elasticity attention maps: additional global vs local examples
Figure 25.Additional attention-map examples on Elasticity, page 2.
52

Hyperbolic Neural Operator
10 4
10 3
10 2
Attention weight (log scale)
Global token (latent 60, d=0.61, RF=0.30)
 Local token (latent 21, d=0.86, RF=0.16)
0.0 0.2 0.4 0.6 0.8 1.0
10 4
10 3
10 2
 Global (d=0.61)
Local (d=0.86)
Global token (latent 48, d=0.66, RF=0.30)
 Local token (latent 7, d=0.86, RF=0.16)
0.0 0.2 0.4 0.6 0.8 1.0
10 4
10 3
10 2
 Global (d=0.66)
Local (d=0.86)
Global token (latent 61, d=0.71, RF=0.31)
 Local token (latent 19, d=0.86, RF=0.22)
0.0 0.2 0.4 0.6 0.8 1.0
10 4
10 3
Global (d=0.71)
Local (d=0.86)
Global token (latent 46, d=0.72, RF=0.26)
 Local token (latent 39, d=0.85, RF=0.19)
0.0 0.2 0.4 0.6 0.8 1.0
Physical distance to barycenter x b 2
10 4
10 3
10 2
Global (d=0.72)
Local (d=0.85)
Mean attention weight
Elasticity attention maps: additional global vs local examples
Figure 26.Additional attention-map examples on Elasticity, page 3.
53

Hyperbolic Neural Operator
10 4
10 3
10 2
Attention weight (log scale)
Global token (latent 53, d=0.73, RF=0.32)
 Local token (latent 38, d=0.85, RF=0.16)
0.0 0.2 0.4 0.6 0.8 1.0
10 4
10 3
10 2
 Global (d=0.73)
Local (d=0.85)
Global token (latent 35, d=0.73, RF=0.25)
 Local token (latent 34, d=0.85, RF=0.16)
0.0 0.2 0.4 0.6 0.8 1.0
10 4
10 3
10 2
 Global (d=0.73)
Local (d=0.85)
Global token (latent 24, d=0.74, RF=0.24)
 Local token (latent 12, d=0.85, RF=0.15)
0.0 0.2 0.4 0.6 0.8 1.0
10 4
10 3
10 2
 Global (d=0.74)
Local (d=0.85)
Global token (latent 17, d=0.75, RF=0.25)
 Local token (latent 32, d=0.85, RF=0.25)
0.0 0.2 0.4 0.6 0.8 1.0
Physical distance to barycenter x b 2
10 4
10 3
10 2
Global (d=0.75)
Local (d=0.85)
Mean attention weight
Elasticity attention maps: additional global vs local examples
Figure 27.Additional attention-map examples on Elasticity, page 4.
54

Hyperbolic Neural Operator
10 4
10 3
10 2
Attention weight (log scale)
Global token (latent 58, d=0.76, RF=0.33)
 Local token (latent 49, d=0.85, RF=0.32)
0.2 0.4 0.6 0.8
10 3
Global (d=0.76)
Local (d=0.85)
Global token (latent 63, d=0.76, RF=0.32)
 Local token (latent 1, d=0.85, RF=0.16)
0.0 0.2 0.4 0.6 0.8 1.0
10 4
10 3
10 2
Global (d=0.76)
Local (d=0.85)
Global token (latent 8, d=0.77, RF=0.21)
 Local token (latent 52, d=0.85, RF=0.19)
0.0 0.2 0.4 0.6 0.8 1.0
10 4
10 3
10 2
 Global (d=0.77)
Local (d=0.85)
Global token (latent 36, d=0.77, RF=0.26)
 Local token (latent 10, d=0.84, RF=0.16)
0.0 0.2 0.4 0.6 0.8 1.0
Physical distance to barycenter x b 2
10 4
10 3
10 2
Global (d=0.77)
Local (d=0.84)
Mean attention weight
Elasticity attention maps: additional global vs local examples
Figure 28.Additional attention-map examples on Elasticity, page 5.
55

Hyperbolic Neural Operator
10 4
10 3
10 2
Attention weight (log scale)
Global token (latent 51, d=0.77, RF=0.32)
 Local token (latent 33, d=0.84, RF=0.19)
0.0 0.2 0.4 0.6 0.8 1.0
10 4
10 3
10 2
 Global (d=0.77)
Local (d=0.84)
Global token (latent 13, d=0.78, RF=0.22)
 Local token (latent 6, d=0.84, RF=0.23)
0.0 0.2 0.4 0.6 0.8 1.0
10 4
10 3
10 2
 Global (d=0.78)
Local (d=0.84)
Global token (latent 47, d=0.79, RF=0.19)
 Local token (latent 3, d=0.84, RF=0.17)
0.0 0.2 0.4 0.6 0.8 1.0
10 4
10 3
10 2
 Global (d=0.79)
Local (d=0.84)
Global token (latent 59, d=0.79, RF=0.20)
 Local token (latent 15, d=0.83, RF=0.15)
0.0 0.2 0.4 0.6 0.8 1.0
Physical distance to barycenter x b 2
10 4
10 3
10 2
Global (d=0.79)
Local (d=0.83)
Mean attention weight
Elasticity attention maps: additional global vs local examples
Figure 29.Additional attention-map examples on Elasticity, page 6.
56

Hyperbolic Neural Operator
10 4
10 3
10 2
Attention weight (log scale)
Global token (latent 54, d=0.79, RF=0.24)
 Local token (latent 16, d=0.83, RF=0.19)
0.0 0.2 0.4 0.6 0.8 1.0
10 4
10 3
10 2
Global (d=0.79)
Local (d=0.83)
Global token (latent 25, d=0.80, RF=0.22)
 Local token (latent 62, d=0.83, RF=0.16)
0.0 0.2 0.4 0.6 0.8 1.0
10 4
10 3
10 2
 Global (d=0.80)
Local (d=0.83)
Global token (latent 43, d=0.81, RF=0.20)
 Local token (latent 29, d=0.83, RF=0.20)
0.0 0.2 0.4 0.6 0.8 1.0
10 4
10 3
10 2
Global (d=0.81)
Local (d=0.83)
Global token (latent 26, d=0.81, RF=0.27)
 Local token (latent 2, d=0.83, RF=0.25)
0.0 0.2 0.4 0.6 0.8 1.0
Physical distance to barycenter x b 2
10 4
10 3
Global (d=0.81)
Local (d=0.83)
Mean attention weight
Elasticity attention maps: additional global vs local examples
Figure 30.Additional attention-map examples on Elasticity, page 7.
57

Hyperbolic Neural Operator
10 4
10 3
10 2
Attention weight (log scale)
Global token (latent 28, d=0.81, RF=0.18)
 Local token (latent 23, d=0.83, RF=0.22)
0.0 0.2 0.4 0.6 0.8 1.0
10 4
10 3
10 2
 Global (d=0.81)
Local (d=0.83)
Global token (latent 0, d=0.81, RF=0.17)
 Local token (latent 42, d=0.82, RF=0.21)
0.0 0.2 0.4 0.6 0.8 1.0
10 4
10 3
10 2
Global (d=0.81)
Local (d=0.82)
Global token (latent 45, d=0.82, RF=0.23)
 Local token (latent 9, d=0.82, RF=0.29)
0.0 0.2 0.4 0.6 0.8 1.0
10 4
10 3
Global (d=0.82)
Local (d=0.82)
Global token (latent 44, d=0.82, RF=0.29)
 Local token (latent 37, d=0.82, RF=0.23)
0.0 0.2 0.4 0.6 0.8
Physical distance to barycenter x b 2
10 4
10 3
10 2
Global (d=0.82)
Local (d=0.82)
Mean attention weight
Elasticity attention maps: additional global vs local examples
Figure 31.Additional attention-map examples on Elasticity, page 8.
58

Hyperbolic Neural Operator
K.3. Additional Mechanism Verification Statistics
We report robust correlations for Fig. 9.Definitions.For a latent query token i with attention weights {αij}N
j=1 over input
nodes at physical coordinates {xj}N
j=1, we define the normalized entropy as Hi/logN , where Hi :=− PN
j=1 αij log(αij).
We define the barycenter bi :=PN
j=1 αijxj and the span (receptive field) RFi :=PN
j=1 αij∥xj −b i∥2, where ∥ · ∥2 is the
Euclidean distance in the physical domain. For each test sample, we compute Spearman correlations between radius r and
either entropy H or physical span RF, flattening over all heads and latent tokens. We also compute head-wise Spearman
correlations within each head (over latent tokens) and report the head-wise median medh ρh per sample. Tab. 17 reports
mean±std over 100 test samples, with 95% bootstrap confidence intervals for the mean.
Table 17.Robust correlations for mechanism verification on Elasticity. Values are mean±std over 100 test samples; bracketed intervals
are 95% bootstrap confidence intervals for the mean.
SCOREρ(r, H)ρ(r,RF)
medh ρh
(r, H)
medh ρh
(r,RF)
HYPERBOLIC
−0.614±0.005
[−0.614,−0.613]
−0.564±0.013
[−0.567,−0.562]
−0.784±0.018
[−0.787,−0.780]
−0.770±0.038
[−0.777,−0.762]
DOT-PRODUCT
(NORMALIZED)
−0.191±0.004
[−0.192,−0.191]
−0.073±0.018
[−0.076,−0.069]
−0.373±0.009
[−0.375,−0.371]
−0.350±0.044
[−0.358,−0.341]
-1
0ρh(r, H)
3 5 8 7 1 4 2 6
head (global → local)
-1
0ρh(r, RF)
1.25 1.31 1.55 1.69 1.94 2.10 2.12 2.51
median radius ̃rh
Hyperbolic Control (cosine)
Figure 32.Head-wise robustness for Fig. 9. Heads are ordered by median query radius (global →local). Markers show the mean and
vertical bars show ±std of per-sample Spearman correlations within each head (over latent tokens); faint points show individual test
samples. Hyperbolic attention yields consistently negative correlations, while the cosine-control is weaker and more heterogeneous.
59