Taylor-Based Hypergraph Inference with SINDy

• The paper proposes a model-free method that fuses a local Taylor expansion with the SINDy algorithm to reconstruct directed hypergraphs from time-series data.
• It identifies higher-order interactions by mapping nonzero Taylor coefficients of mixed monomials to corresponding hyperedges.
• The approach is extended with Bayesian techniques (Bayes-THIS) to quantify uncertainty and address structural non-identifiability challenges.

Searching arXiv for THIS, Bayes-THIS, and related SINDy hypergraph inference papers. Taylor-based Hypergraph Inference using SINDy (THIS) is a model-free hypergraph inference method that reconstructs higher-order interaction structure from time-series data by combining a local Taylor expansion of node dynamics with sparse identification of nonlinear dynamics (SINDy). In this framework, nonzero Taylor coefficients associated with mixed monomials are interpreted as evidence of pairwise or higher-order interactions, allowing reconstruction of directed hypergraphs and simplicial complexes from dynamical observations without requiring prior knowledge of the node dynamics or coupling functions (Delabays et al., 2024). Subsequent work has positioned THIS as both an inference method in its own right and an identification stage within data-driven control pipelines for diffusive hypergraph dynamical systems, while also clarifying its statistical limitations and motivating Bayesian extensions such as Bayes-THIS (Delabays et al., 10 Nov 2025, Tang et al., 5 May 2026).

1. Conceptual basis and problem formulation

THIS addresses the problem of inferring nonpairwise interactions from observed dynamics. The motivating setting is a coupled ordinary differential equation system on $n$ nodes, where node $i$ evolves according to

$$\dot{x}_i = F_i(\mathbf{x}) = f_i(x_i) + \sum_{j=1}^n a_{ij}^{(2)} h^{(2)}(x_i,x_j) + \sum_{j,k=1}^n a_{ijk}^{(3)} h^{(3)}(x_i,x_j,x_k)+\cdots$$

or, in an equivalent notation used later,

$$\dot{x}_i = F_i(\mathbf{x}) = f_i(x_i) + \sum_{j=1}^n a_{ij}^{(2)} g^{(2)}(x_i,x_j) + \sum_{j,k=1}^n a_{ijk}^{(3)} g^{(3)}(x_i,x_j,x_k) + \cdots .$$

Here, $a^{(2)}$ encodes pairwise edges and $a^{(3)}$ encodes triadic or hyperedges (Delabays et al., 2024, Tang et al., 5 May 2026).

The method is motivated by the observation that many real-world systems cannot be adequately captured by ordinary pairwise graphs because the state of a node may depend simultaneously on multiple other nodes. The control-oriented formulation in later work emphasizes directed hypergraphs, in which a $k$-hyperedge is represented by a tuple

$e=(i_1,i_2,\dots,i_k),$

with node $i_1$ as the head and nodes $i_2,\dots,i_k$ as the tails, so that the tails jointly influence the head. The associated structural data are encoded in higher-order adjacency tensors $i$0, where

$i$1

indicates that the directed $i$2-hyperedge is present (Delabays et al., 10 Nov 2025).

Within this formalism, THIS reconstructs the incoming hyperedges to each node independently. For each target node $i$3, it solves a sparse regression problem over a library of candidate nonlinear terms evaluated on observed trajectories. In the original formulation, the output is binary—edge present or absent—rather than a weighted hyperedge estimate (Delabays et al., 2024). Later work emphasizes that the inferred sparse coefficients also encode interaction strengths in addition to structural support (Delabays et al., 10 Nov 2025). This suggests that the operational interpretation of THIS depends on the downstream task: in topology reconstruction, support recovery is primary, whereas in control-oriented settings the coefficient magnitudes may also be used as effective local coupling strengths.

2. Taylor expansion and the SINDy regression architecture

The defining step of THIS is the replacement of unknown interaction functions by a local Taylor approximation around a reference point. In the original reconstruction framework, the dynamics of node $i$4 are approximated around an arbitrary point $i$5 as

$i$6

with

$i$7

A later formulation writes the same idea around a reference point $i$8 as

$i$9

where $$\dot{x}_i = F_i(\mathbf{x}) = f_i(x_i) + \sum_{j=1}^n a_{ij}^{(2)} h^{(2)}(x_i,x_j) + \sum_{j,k=1}^n a_{ijk}^{(3)} h^{(3)}(x_i,x_j,x_k)+\cdots$$0 (Delabays et al., 2024, Tang et al., 5 May 2026).

The central interpretive rule is that a monomial of degree $$\dot{x}_i = F_i(\mathbf{x}) = f_i(x_i) + \sum_{j=1}^n a_{ij}^{(2)} h^{(2)}(x_i,x_j) + \sum_{j,k=1}^n a_{ijk}^{(3)} h^{(3)}(x_i,x_j,x_k)+\cdots$$1 in the Taylor expansion corresponds to an interaction of order $$\dot{x}_i = F_i(\mathbf{x}) = f_i(x_i) + \sum_{j=1}^n a_{ij}^{(2)} h^{(2)}(x_i,x_j) + \sum_{j,k=1}^n a_{ijk}^{(3)} h^{(3)}(x_i,x_j,x_k)+\cdots$$2. Linear terms correspond to pairwise structure, quadratic terms correspond to triadic or 3-body interactions, and cubic terms correspond to 4-body interactions (Delabays et al., 2024). In the Bayesian exposition, the interpretation is stated in derivative form: linear terms $$\dot{x}_i = F_i(\mathbf{x}) = f_i(x_i) + \sum_{j=1}^n a_{ij}^{(2)} h^{(2)}(x_i,x_j) + \sum_{j,k=1}^n a_{ijk}^{(3)} h^{(3)}(x_i,x_j,x_k)+\cdots$$3 correspond to pairwise effects, while quadratic terms $$\dot{x}_i = F_i(\mathbf{x}) = f_i(x_i) + \sum_{j=1}^n a_{ij}^{(2)} h^{(2)}(x_i,x_j) + \sum_{j,k=1}^n a_{ijk}^{(3)} h^{(3)}(x_i,x_j,x_k)+\cdots$$4 correspond to triadic or higher-order effects (Tang et al., 5 May 2026). In particular, if $$\dot{x}_i = F_i(\mathbf{x}) = f_i(x_i) + \sum_{j=1}^n a_{ij}^{(2)} h^{(2)}(x_i,x_j) + \sum_{j,k=1}^n a_{ijk}^{(3)} h^{(3)}(x_i,x_j,x_k)+\cdots$$5 are distinct and $$\dot{x}_i = F_i(\mathbf{x}) = f_i(x_i) + \sum_{j=1}^n a_{ij}^{(2)} h^{(2)}(x_i,x_j) + \sum_{j,k=1}^n a_{ijk}^{(3)} h^{(3)}(x_i,x_j,x_k)+\cdots$$6, this indicates a genuine triadic interaction involving nodes $$\dot{x}_i = F_i(\mathbf{x}) = f_i(x_i) + \sum_{j=1}^n a_{ij}^{(2)} h^{(2)}(x_i,x_j) + \sum_{j,k=1}^n a_{ijk}^{(3)} h^{(3)}(x_i,x_j,x_k)+\cdots$$7, assuming truncation to at most third-order interactions (Tang et al., 5 May 2026).

SINDy supplies the estimation machinery. The original paper writes

$$\dot{x}_i = F_i(\mathbf{x}) = f_i(x_i) + \sum_{j=1}^n a_{ij}^{(2)} h^{(2)}(x_i,x_j) + \sum_{j,k=1}^n a_{ijk}^{(3)} h^{(3)}(x_i,x_j,x_k)+\cdots$$8

where $$\dot{x}_i = F_i(\mathbf{x}) = f_i(x_i) + \sum_{j=1}^n a_{ij}^{(2)} h^{(2)}(x_i,x_j) + \sum_{j,k=1}^n a_{ijk}^{(3)} h^{(3)}(x_i,x_j,x_k)+\cdots$$9 is a library or design matrix of nonlinear monomials evaluated on the data and $$\dot{x}_i = F_i(\mathbf{x}) = f_i(x_i) + \sum_{j=1}^n a_{ij}^{(2)} g^{(2)}(x_i,x_j) + \sum_{j,k=1}^n a_{ijk}^{(3)} g^{(3)}(x_i,x_j,x_k) + \cdots .$$0 is a sparse coefficient vector (Delabays et al., 2024). Closely related notation used in later work is

$$\dot{x}_i = F_i(\mathbf{x}) = f_i(x_i) + \sum_{j=1}^n a_{ij}^{(2)} g^{(2)}(x_i,x_j) + \sum_{j,k=1}^n a_{ijk}^{(3)} g^{(3)}(x_i,x_j,x_k) + \cdots .$$1

or, for a single node,

$$\dot{x}_i = F_i(\mathbf{x}) = f_i(x_i) + \sum_{j=1}^n a_{ij}^{(2)} g^{(2)}(x_i,x_j) + \sum_{j,k=1}^n a_{ijk}^{(3)} g^{(3)}(x_i,x_j,x_k) + \cdots .$$2

with $$\dot{x}_i = F_i(\mathbf{x}) = f_i(x_i) + \sum_{j=1}^n a_{ij}^{(2)} g^{(2)}(x_i,x_j) + \sum_{j,k=1}^n a_{ijk}^{(3)} g^{(3)}(x_i,x_j,x_k) + \cdots .$$3 the library matrix and $$\dot{x}_i = F_i(\mathbf{x}) = f_i(x_i) + \sum_{j=1}^n a_{ij}^{(2)} g^{(2)}(x_i,x_j) + \sum_{j,k=1}^n a_{ijk}^{(3)} g^{(3)}(x_i,x_j,x_k) + \cdots .$$4 sparse (Delabays et al., 10 Nov 2025, Tang et al., 5 May 2026).

The candidate library is a monomial basis. One explicit description is

$$\dot{x}_i = F_i(\mathbf{x}) = f_i(x_i) + \sum_{j=1}^n a_{ij}^{(2)} g^{(2)}(x_i,x_j) + \sum_{j,k=1}^n a_{ijk}^{(3)} g^{(3)}(x_i,x_j,x_k) + \cdots .$$5

while another writes

$$\dot{x}_i = F_i(\mathbf{x}) = f_i(x_i) + \sum_{j=1}^n a_{ij}^{(2)} g^{(2)}(x_i,x_j) + \sum_{j,k=1}^n a_{ijk}^{(3)} g^{(3)}(x_i,x_j,x_k) + \cdots .$$6

Each coefficient multiplies one monomial term in the Taylor basis, and nonzero coefficients are mapped to pairwise or higher-order hyperedges (Delabays et al., 2024, Tang et al., 5 May 2026).

The sparse solver in the original THIS uses SINDy with sequential thresholded least squares. The Bayesian discussion notes that Delabays et al.’s original THIS uses sparse regression methods such as LASSO or sequential thresholded least squares, and that sequential thresholded least squares employs a global hard-thresholding parameter $$\dot{x}_i = F_i(\mathbf{x}) = f_i(x_i) + \sum_{j=1}^n a_{ij}^{(2)} g^{(2)}(x_i,x_j) + \sum_{j,k=1}^n a_{ijk}^{(3)} g^{(3)}(x_i,x_j,x_k) + \cdots .$$7 (Delabays et al., 2024, Tang et al., 5 May 2026). This places THIS within the broader SINDy lineage, where a large candidate library is pruned to a parsimonious active set. The more general SINDy literature also provides related extensions, such as implicit-SINDy for rational nonlinearities, which likewise interprets active mixed monomials as interaction structure (Mangan et al., 2016).

3. Hypergraph representation, inference targets, and identifiability conditions

THIS is designed to reconstruct directed hypergraphs, simplicial complexes, and general higher-order hypergraphs with interactions up to a chosen maximal order (Delabays et al., 2024). In the directed hypergraph convention used in the control paper, the first tensor index is the head and the remaining indices are the tails. Thus, a nonzero coefficient associated with a candidate term involving nodes $$\dot{x}_i = F_i(\mathbf{x}) = f_i(x_i) + \sum_{j=1}^n a_{ij}^{(2)} g^{(2)}(x_i,x_j) + \sum_{j,k=1}^n a_{ijk}^{(3)} g^{(3)}(x_i,x_j,x_k) + \cdots .$$8 corresponds to a directed $$\dot{x}_i = F_i(\mathbf{x}) = f_i(x_i) + \sum_{j=1}^n a_{ij}^{(2)} g^{(2)}(x_i,x_j) + \sum_{j,k=1}^n a_{ijk}^{(3)} g^{(3)}(x_i,x_j,x_k) + \cdots .$$9-hyperedge from $a^{(2)}$0 to $a^{(2)}$1 (Delabays et al., 10 Nov 2025).

The distinction between pairwise and nonpairwise interactions is made through irreducible mixed terms in the Taylor representation. The original reconstruction paper notes that

$a^{(2)}$2

whereas

$a^{(2)}$3

is generally nonzero. On this basis, a nonzero mixed partial derivative is evidence of a nonpairwise interaction (Delabays et al., 2024). However, the same paper also states a major ambiguity: if the triadic interaction function $a^{(2)}$4 can be written as a linear combination of pairwise interaction functions, then there is no formal difference between a 2-simplex and a closed triangle of links (Delabays et al., 2024). This is a fundamental identifiability condition: decomposable higher-order interactions are not distinguishable from pairwise networks within the dynamical observations alone.

A second identifiability condition concerns the sampled region of state space. The original paper states that if data only come from a small neighborhood in state space, nonlinear interactions may appear effectively linear, and hypergraph inference becomes impossible (Delabays et al., 2024). The later control paper makes the same point operationally: the explored region should be large enough to reveal nonlinearity beyond a purely linear neighborhood, but not so large that the Taylor approximation breaks down (Delabays et al., 10 Nov 2025).

The original paper further notes that if the system is known to be a simplicial complex, then the downward inclusion condition can be used to infer lower-order interactions automatically. It also states that the same idea works approximately for anti-simplicial complexes and random sparse hypergraphs where overlap between lower- and higher-order interactions is limited (Delabays et al., 2024). This suggests that prior structural assumptions can materially affect interpretation of the inferred support, especially when lower-order terms may arise from higher-order mechanisms.

4. Data requirements and practical workflow

The basic workflow of THIS is consistent across its original and subsequent formulations. The method begins with time-series observations $a^{(2)}$5, forms or estimates derivatives, constructs a Taylor-appropriate monomial library, solves a sparse regression problem for each node, and maps active monomials to hyperedges (Delabays et al., 2024, Delabays et al., 10 Nov 2025). In the original benchmarks, exact derivatives were computed from the ODE in some experiments, while numerical differentiation was used in others; the paper states that finite differences can be used in practice and that robust SINDy variants like Ensemble-SINDy and Weak-SINDy could improve performance (Delabays et al., 2024).

The control-oriented deployment of THIS makes the data regime particularly explicit. There, the system is initialized near a locally stable equilibrium $a^{(2)}$6, allowed to evolve without control over a time window $a^{(2)}$7, and the resulting noisy trajectories are used as data input to THIS (Delabays et al., 10 Nov 2025). The noise is not merely a nuisance. The paper states that the noisy system passively explores a neighborhood of the equilibrium, and that this exploration is helpful because it reveals local nonlinear structure without requiring active system perturbation (Delabays et al., 10 Nov 2025).

The original reconstruction paper likewise emphasizes the importance of the sampling region. Data were sampled from a box around the origin in the Kuramoto benchmark, and the authors recommend scanning a range of box sizes and looking for a plateau of stable performance because too small a region yields nearly linear dynamics whereas too large a region degrades the Taylor approximation (Delabays et al., 2024). A plausible implication is that THIS is best viewed as a local structural inference method whose success depends on choosing an observational neighborhood that balances expressiveness against truncation error.

Scalability is limited by the combinatorial growth of the monomial library. The original paper states that the algorithm currently scales roughly like about $a^{(2)}$8 in the tested third-order setting, though parallelization across nodes can improve practical scaling. It also recommends pre-filtering node pairs by correlation to reduce library size in large systems (Delabays et al., 2024). For incomplete observations, the same paper states that THIS should be combined with data assimilation methods such as ensemble Kalman filtering (Delabays et al., 2024).

5. Empirical benchmarks and applications

The original empirical validation of THIS focused on synthetic dynamical systems with known higher-order structure. In a higher-order Kuramoto benchmark,

$a^{(2)}$9

the method was tested first on a 7-node hypergraph and then on larger systems including a 100-node simplicial complex, a 100-node random hypergraph, and a 300-node random simplicial complex with filtering (Delabays et al., 2024). Performance was measured by ROC curves, and the paper reports that with just 10 data points, THIS already reconstructs the hypergraph well. It also states that ARNI is much weaker here unless the function basis is chosen very carefully, while for pairwise-only inference ARNI can match THIS if the right basis is chosen (Delabays et al., 2024).

A second benchmark used coupled Lorenz systems with nonpairwise coupling: $a^{(3)}$0 with

$a^{(3)}$1

For 5-node random hypergraphs, the paper reports ROC curves showing over 80% TPR with less than 20% FPR in many cases, and states that performance improves with more data, better box-size tuning, or exact derivatives (Delabays et al., 2024).

The most prominent real-data application in the original paper is resting-state EEG. The dataset comprises 109 subjects with 2 recordings each, for a total of 218 time series. Starting from 64 sensors, the signals were coarse-grained to 7 brain regions, low-pass filtered to below about 5 Hz, standardized to unit variance, restricted to the 1000 data points closest to the median, and differentiated by finite differences (Delabays et al., 2024). THIS was used to infer effective connectivity among the seven brain regions, and the paper reports that nonpairwise interactions account for more than 60% of the EEG dynamics. Across thresholds, pairwise terms explain about 35%, third-order terms about 45%, and fourth-order terms about 20% (Delabays et al., 2024). The most frequent 3-edges and 4-edges point toward the region corresponding roughly to the prefrontal cortex (Delabays et al., 2024).

The following table summarizes the benchmark and application settings explicitly described in the literature.

Setting	System	Reported finding
Synthetic benchmark	Higher-order Kuramoto model	With just 10 data points, THIS already reconstructs the hypergraph well (Delabays et al., 2024)
Synthetic benchmark	Coupled Lorenz systems	ROC curves show over 80% TPR with less than 20% FPR in many cases (Delabays et al., 2024)
Empirical application	Resting-state EEG	Nonpairwise interactions account for more than 60% of the EEG dynamics (Delabays et al., 2024)

These results established THIS as a practical bridge between local nonlinear approximation, sparse system identification, and higher-order topology inference (Delabays et al., 2024).

6. Control-oriented use in diffusive hypergraph dynamical systems

A later development embeds THIS into a closed-loop control framework for diffusive hypergraph dynamical systems (Delabays et al., 10 Nov 2025). The paper considers systems operating near a stable equilibrium but subject to noise, with dynamics of the form

$a^{(3)}$2

with unknown interaction functions $a^{(3)}$3, where only time-series data $a^{(3)}$4 are available (Delabays et al., 10 Nov 2025). In this setting, the controller cannot be designed from first principles because the true hypergraph is unknown; instead, control must be built on top of an inferred model.

The proposed workflow is: first observe noisy trajectories near a stable equilibrium; then infer which higher-order couplings are present using sparse regression on a Taylor-expanded candidate library; and finally use the inferred structure to choose a minimal set of control nodes and apply a simple proportional controller to suppress noise-driven departures from the desired equilibrium (Delabays et al., 10 Nov 2025). The control design is based on the notion of a leaf node, defined as a node that is not the head of any edge. The minimal controllable set is

$a^{(3)}$5

The rationale given is that if a node is not the head of any hyperedge, it cannot be reached through directed influence from other nodes, so it must be actuated directly; once all leaf nodes are controlled, inputs can propagate through directed paths to the rest of the system (Delabays et al., 10 Nov 2025).

The controller is a parsimonious droop or proportional state-feedback law applied only on nodes in $a^{(3)}$6: $a^{(3)}$7 with $a^{(3)}$8 chosen to preserve stability of the desired equilibrium $a^{(3)}$9 (Delabays et al., 10 Nov 2025). The stated control objective is to design a parsimonious state-feedback controller $k$0 so that the system converges to the closest equilibrium point $k$1 (Delabays et al., 10 Nov 2025).

Validation is carried out on a third-order Kuramoto model on a random 3-hypergraph: $k$2 where $k$3 is the phase of oscillator $k$4, $k$5 its natural frequency, and $k$6 and $k$7 encode pairwise and 3-body couplings (Delabays et al., 10 Nov 2025). The reported example uses 10 oscillators and noisy trajectories near the synchronous state as data. The inferred topology is imperfect: the true positive rate for 3-edges is only 55%, and some spurious 2-edges are inferred even though the true system had none (Delabays et al., 10 Nov 2025). Nonetheless, THIS correctly identifies the leaf nodes, after which the proportional control law with $k$8 is applied only to those leaf nodes starting at $k$9, damping fluctuations and driving the whole network toward synchrony (Delabays et al., 10 Nov 2025).

The principal conclusion of this control study is that exact recovery of every hyperedge is not necessary for effective control. What matters is identification of the critical structural information needed for actuation and reachability, particularly the leaf nodes (Delabays et al., 10 Nov 2025). This shifts the interpretation of THIS from a pure topology estimator toward a control-oriented inference tool.

7. Bayesian extension and fundamental limitations

Bayes-THIS extends THIS by replacing fixed-threshold sparse regression with sparse Bayesian regression using automatic relevance determination (ARD) (Tang et al., 5 May 2026). The representation remains the same: hypergraph structure is inferred by identifying sparse Taylor coefficients associated with pairwise and higher-order interactions. The regression model for one node is

$e=(i_1,i_2,\dots,i_k),$0

with likelihood

$e=(i_1,i_2,\dots,i_k),$1

Each coefficient receives an independent Gaussian ARD prior,

$e=(i_1,i_2,\dots,i_k),$2

with Gamma hyperprior

$e=(i_1,i_2,\dots,i_k),$3

and a minimally informative choice $e=(i_1,i_2,\dots,i_k),$4 (Tang et al., 5 May 2026).

Evidence maximization yields a Gaussian posterior

$e=(i_1,i_2,\dots,i_k),$5

with

$e=(i_1,i_2,\dots,i_k),$6

where $e=(i_1,i_2,\dots,i_k),$7 and $e=(i_1,i_2,\dots,i_k),$8 (Tang et al., 5 May 2026). Compared with sequential thresholded least squares, this yields term-wise shrinkage, residual variance modeling, and posterior uncertainty (Tang et al., 5 May 2026).

Bayes-THIS also introduces an uncertainty-aware workflow. Before full higher-order inference, the library is partitioned as

$e=(i_1,i_2,\dots,i_k),$9

where $i_1$0 contains pairwise monomials and $i_1$1 triadic monomials. A posterior predictive check (PPC) is then used to assess whether the data contain sufficient higher-order signal to support inference beyond a pairwise model (Tang et al., 5 May 2026). If the PPC is favorable, the full model is fit and coefficients are pruned using credible intervals: a term is retained if

$i_1$2

where $i_1$3 is the conditional posterior variance (Tang et al., 5 May 2026). The paper characterizes this as uncertainty-aware and robust in practice, though not fully Bayesian in the strictest sense because it conditions on point estimates for other coefficients and hyperparameters (Tang et al., 5 May 2026).

The most consequential result of the Bayesian paper is not merely algorithmic but diagnostic. It identifies a structural non-identifiability in the Taylor-based inference framework: higher-order interactions can generate lower-order Taylor terms, producing spurious lower-order edges that are indistinguishable from genuine lower-order interactions (Tang et al., 5 May 2026). The paper emphasizes a particularly bad regime of low cross-order degree correlation, in which higher-order interactions occur on nodes that do not already have lower-order connections. In that case, the downward contribution from triadic interactions lands on absent pairwise edges, making them look like real pairwise couplings (Tang et al., 5 May 2026).

This limitation had already appeared empirically in the original reconstruction paper, which noted that if a $i_1$4-edge exists, THIS may also infer some sub-edges even when they do not exist, especially in non-simplicial hypergraphs (Delabays et al., 2024). The control paper provides a concrete instance: in the 10-oscillator Kuramoto example, some spurious 2-edges were inferred even though the true system had none (Delabays et al., 10 Nov 2025). Bayes-THIS sharpens this observation into a representation-level impossibility result: no threshold choice, credible interval, or Bayesian regularization can resolve the ambiguity within the Taylor-based framework alone (Tang et al., 5 May 2026).

Taken together, these studies define the present understanding of THIS. It is a Taylor-expansion-based hypergraph inference method that translates dynamical observations into sparse regression over monomial terms, thereby recovering candidate higher-order structure from data (Delabays et al., 2024). It can function effectively as a control-oriented identification engine even when edge-level recovery is imperfect (Delabays et al., 10 Nov 2025). At the same time, its inference capabilities are constrained by local sampling, algebraic decomposability, and a fundamental structural non-identifiability whereby higher-order interactions can masquerade as lower-order ones (Tang et al., 5 May 2026).