Hyperbolic Neural Operator Jieyuan Pei, Zhuoxuan Li, Wei Li, Haobo Zhang, Jiawei Jiang, Jianwei Zheng ICML 2026 · PMLR 306 Source PDF: https://guobapei.github.io/Hyperbolic-Neural-Operator/assets/paper.pdf Text extracted from the author-provided PDF. Consult the PDF for equations, figures and table layout. --- PDF page 1 --- 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 . 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 --- PDF page 2 --- 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 --- PDF page 3 --- 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 --- PDF page 4 --- 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 --- PDF page 5 --- 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 --- PDF page 6 --- 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 --- PDF page 7 --- 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 --- PDF page 8 --- 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 --- PDF page 9 --- 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 --- PDF page 10 --- 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 --- PDF page 11 --- 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 --- PDF page 12 --- 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 --- PDF page 13 --- 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 --- PDF page 14 --- 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 --- PDF page 15 --- 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 --- PDF page 16 --- 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 --- PDF page 17 --- 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 --- PDF page 18 --- 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 --- PDF page 19 --- 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 --- PDF page 20 --- 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 --- PDF page 21 --- 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 --- PDF page 22 --- 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 --- PDF page 23 --- 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 --- PDF page 24 --- 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 --- PDF page 25 --- 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 --- PDF page 26 --- 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 --- PDF page 27 --- 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 --- PDF page 28 --- 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 --- PDF page 29 --- 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 --- PDF page 30 --- 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 --- PDF page 31 --- 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∥22 +∥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 --- PDF page 32 --- 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 --- PDF page 33 --- 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 --- PDF page 34 --- 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 --- PDF page 35 --- 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 --- PDF page 36 --- 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 --- PDF page 37 --- 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 --- PDF page 38 --- 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 --- PDF page 39 --- 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 --- PDF page 40 --- 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 --- PDF page 41 --- 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 --- PDF page 42 --- 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 --- PDF page 43 --- 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 --- PDF page 44 --- 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 --- PDF page 45 --- 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 --- PDF page 46 --- 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 --- PDF page 47 --- 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 --- PDF page 48 --- 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 --- PDF page 49 --- 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 --- PDF page 50 --- 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 --- PDF page 51 --- 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 --- PDF page 52 --- 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 --- PDF page 53 --- 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 --- PDF page 54 --- 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 --- PDF page 55 --- 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 --- PDF page 56 --- 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 --- PDF page 57 --- 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 --- PDF page 58 --- 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 --- PDF page 59 --- 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