A simplex path integral and a simplex renormalization group for high-order interactions^*

In recent decades, the studies on phase transition phenomena, especially non-equilibrium ones, in different interacting systems have made substantial progress [1, 2]. This progress should be credited to the development of path integrals [3–5] and renormalization group (RG) [6, 7] theories, which significantly deepen our understanding of system dynamics and provide a precise formulation of scaling and criticality.

However, a theoretical vacancy can be found in the existing path integral and RG approaches if we subdivide interactions into pairwise and high-order categories. As the name suggests, a pairwise interaction only involves a pair of units and does not require the participation of any other unit. These kinds of interactions can be represented by edges in networks [8, 9] and are implicitly used in the derivations of classic path integrals and RG theories (e.g. see [10]). A high-order interaction, in contrast, is the mutual coupling among more than two units [11–13]. While some high-order interactions can be decomposed into a group of pairwise interactions, most high-order interactions are undecomposable and in-equivalent to the direct sum of the pairwise interactions [9, 14]. For instance, the triplet collaboration among three agents is not equivalent to the trivial sum of three pairs of individual collaborations. These undecomposable high-order interactions have intricate effects on system dynamics [13, 14] and cannot be trivially represented by ordinary networks [12]. Although notable efforts have been devoted to studying the optimal characterization of high-order interactions (e.g. using simplicial complexes [15–18] or hypergraphs [19–21]), there are fewer works focusing on the development of specialized path integrals and RG theories for complex systems with high-order interactions [17, 22]. The difficulty in proposing these specialized frameworks arises from the intrinsic dependence on pairwise interactions when researchers constrain the discrete replicas of interacting systems as conventional networks. Once this constraint is relaxed, the intertwined effects of newly included high-order interactions would make the analysis of emerging dynamics highly non-trivial and even lead to novel results that cannot be directly derived by classic theories (e.g. see instances in diffusion and random walks [16, 20, 23], social dynamics [24–26], and neural dynamics [27, 28]).

To suggest a possible direction for developing path integrals and RG theories for high-order interactions, we summarize the accomplishments and limitations of previous works.

To date, most of the progress of the computational implementations of path integrals and RGs has been achieved in real systems with pairwise interactions [10, 29–33]. Although box-covering methods [34–37] may be the most straightforward methods for coarse graining while preserving the organizational properties of a system represented by a network, they depend on the assumption of internal fractal properties and convey no information about the system dynamics. In contrast to box-covering, the spectral method [38] pays more attention to inherent dynamics properties (e.g. random walks) during coarse graining but may suffer from the obstacles of finding macro-units implied by small-world effects [10, 33, 39]. With the aim to resolve the coexisting and correlated scales, a geometric RG has been developed to identify the potential geometric scaling rather than path-distance scaling [33]. Nevertheless, this approach essentially relies on the a priori knowledge of the hidden metric space of networks (e.g. the hyperbolic space [40]). With the pursuit of realizing metric-free coarse graining, whose iterability is not limited by small-world effects, a Laplacian renormalization group (LRG) [10] has been derived from the diffusion on networks, a continuum counterpart of the block-spin transformation [41], to deal with the network topology, dynamics, and geometry encoded by the graph Laplacian operator [42]. This approach has a natural connection with path integrals as well [43]. Although these properties make the LRG a promising foundation for in-depth explorations, this framework suffers from several non-negligible limitations. First, there is no apparent correspondence between the LRG and an appropriate approach for analyzing undecomposable high-order interactions because the Laplacian operator cannot characterize polyadic relations directly [9, 11–14]. Second, the intrinsic dependence of the LRG on the ergodic properties of the system dynamics reflected by network connectivity limits the applicability of this approach to analyzing high-order interactions, which are more sparsely distributed in the system than pairwise interactions and do not necessarily ensure ergodicity.

Fewer pioneering works exist that are devoted to developing RG theories for high-order interactions, among which a notable framework is the real space RG proposed by [17, 22]. This approach is initially defined on the Laplacian operators of the network skeletons of the Apollonian [44] and the pseudo-fractal [45] simplicial complexes to calculate spectral dimensions [22]. Then, it is generalized to the normalized up-Laplacian operator [16, 23] of simplicial complexes to deal with high-order spectra. Despite the effectiveness of this approach, its dependence on the non-trivial manual derivations using the Gaussian model [46] and the specialized application scope of spectral dimensions [17, 22] make it less applicable to real scenarios, where a programmatic implementation is demanded for renormalization of complex systems generated by empirical data automatically.

In summary, while representative works such as the LRG [10] have suggested a promising way to renormalize complex systems, the appropriate path integrals and RG theories for high-order interactions remain clouded by numerous theoretical unknowns and technical difficulties.

To fill these gaps, we generalize classic path integrals and RG theories from dyadic to polyadic relations. To realize such generalizations, we are required to derive the appropriate definitions of the canonical density operators (i.e. the functional over fields) [10, 49, 50] and the coarse graining processes in real and moment spaces [6, 7] for high-order interactions.

Our work suggests a way to define these concepts in simplices, in which simplices can join any subset of units to characterize undecomposable high-order interactions among them (see figure 1(a) for instances and section I in the supplementary materials for definitions). We propose two kinds of high-order propagators as the generalized path integral formulation of the abstract diffusion processes driven by multi-order interactions, which further creates a system mode description tantamount to the moment space. The coarse graining in moment space is implemented by progressively integrating out short-range system modes corresponding to certain orders of interactions. Parallel to this process, we propose a real-space coarse graining procedure to reduce the simplex structures associated with short-range high-order interactions, where we suggest a way to overcome the dependence on ergodicity. The correspondence between the moment and real spaces enables our framework to implement system coarse graining according to both organizational structure and internal dynamics. Taken together, these definitions generate our simplex path integral and simplex renormalization group (SRG) theories, which are applicable to arbitrary orders of interactions.

Zoom In Zoom Out Reset image size

The SRG can be used to uncover the potential scale invariance of high-order interactions, whose non-trivial fixed point may suggest the existence of critical phenomena. In contrast to conventional RGs, the SRG is designed to renormalize p-order interactions (each p-order interaction requires the simultaneous participation of p + 1 units) in the system according to the structure and dynamics associated with q-order interactions ($p\unicode{x2A7D} q$). Therefore, we can compare the renormalization flow of p-order interactions guided by p-order interactions with that guided by arbitrary higher-order interactions. This feature allows us to investigate the impacts of higher-order interactions on lower-order interactions and study the differences between lower-order and higher-order interactions in terms of scaling properties. In a special case with $p = q = 1$ (i.e. pairwise interactions), the two kinds of high-order propagators developed in our work both reduce to the classic counterpart used in the LRG [10], and our proposed SRG can reproduce the outputs of the LRG in the ergodic case.

To comprehensively demonstrate the capability of our approaches to identify the latent characteristics of diverse real complex systems, we validate our framework from the perspectives of multi-order scale-invariance detection, topological invariance discovery, latent organizational structure identification, and information bottleneck optimization, respectively. An efficient code implementation of this framework is provided in [51]. One can also see a detailed explanation about the programmatic implementation and usage of the SRG in section XIII in the supplementary materials.

Let us consider the n-order interactions that manifest as the simultaneous participation of n + 1 related units. The undecomposability of these interactions lies in the fact that the couplings among n + 1 units are required to occur simultaneously and none of them are dispensable. For instance, the triplet collaboration requires the participation of all three agents together. As discussed in [9], these interactions can be represented by the n-faces of the n-simplex, where n + 1 units are placed. For the sake of clarity, we define $\{i\}_{n+1}\subset V$ as an arbitrary set of n + 1 units in which unit i is included, and we denote $S^{\{i\}}_{n+1}$ as the set of all permutations (i.e. rearrangements) on $\{i\}_{n+1}$. The multi-order Laplacian $\mathbf{L}_{M}^{\left(n\right)}$ associated with n-simplices is defined as

$$\begin{align} \left[\mathbf{L}_{M}^{\left(n\right)}\right]_{ij} = n \delta_{ij}\mathbf{D}_{i}^{\left(n\right)}-\mathbf{A}_{ij}^{\left(n\right)}.\end{align} \tag{ 1 }$$

In equation (1), the term $\mathbf{D}_{i}^{\left(n\right)}$ counts the number of n-simplices that contains unit i

$$\begin{align} \mathbf{D}_{i}^{\left(n\right)} = \sum_{\left\{i\right\}_{n+1}\subset V}\prod_{\omega\in S^{\left\{i\right\}}_{n+1}}\mathbf{M}_{\omega}.\end{align} \tag{ 2 }$$

Meanwhile, the term $\mathbf{A}_{ij}^{\left(n\right)}$ counts the number of n-simplices that contains units i and j

$$\begin{align} \mathbf{A}_{ij}^{\left(n\right)} = \left(1-\delta_{ij}\right)\sum_{\left\{i,j\right\}_{n+1}\subset V}\prod_{\omega\in S^{\left\{i,j\right\}}_{n+1}}\mathbf{M}_{\omega},\end{align} \tag{ 3 }$$

where $\{i,j\}_{n+1}\subset V$ is an arbitrary set of n + 1 units in which units i and j are included. The general definition of $\mathbf{M}_{\omega}$ used in equations (2) and (3) is given as a product of a series of matrices for any permutation $\omega = \left(k_{1}\ldots k_{n}\right)$

$$\begin{align} \mathbf{M}_{\omega} = \mathbf{A}_{k_{1}k_{2}}\cdots\mathbf{A}_{k_{n-1}k_{n}},\end{align} \tag{ 4 }$$

where $\mathbf{A} = \mathbf{A}^{\left(1\right)}$ is the classic adjacency matrix of units. It can be seen that $\prod_{\omega\in S^{\{i,j\}}_{n+1}}\mathbf{M}_{\omega}$ serves as an indicator function, which equals 1 if the units in $\{i,j\}_{n+1}$ form an n-simplex and equals 0 otherwise.

See section II in the supplementary materials for the full derivations of the multi-order Laplacian. In equation (2), vector $\mathbf{D}^{\left(n\right)}$ is similar to the degree vector, while matrix $\mathbf{A}^{\left(n\right)}$ can be understood as the adjacency matrix. The coefficient n in equation (2) arises from the fact that $\sum_{j}\mathbf{A}_{ij}^{\left(n\right)} = n\mathbf{D}^{\left(n\right)}_{i}$ (i.e. each simplex is repeatedly counted n times). In a special case where n = 1, equation (2) is equivalent to the classic Laplacian for pairwise interactions. In the more general cases where n > 1, high-order interactions can be characterized by equation (2). See figure 1(d) for illustrations of the multi-order Laplacian. Although the definitions of equations (2)–(4) are different from the approach proposed by [9], they are mathematically consistent with [9].

Then, we turn to another type of n-order interactions, which manifest as the sequential actions of n + 1 related units after they all participate in (i.e. sequential participation). The undecomposability of these interactions arises from the fact that the mutual coupling effects emerge only after n + 1 units have been engaged with the interplay progressively. For example, a pipeline-like collaboration may require all agents to communicate with each other and behave in order. Each previous action is globally known to all agents and affects all subsequent actions. This property can be characterized by the paths (i.e. a sequence of 1-simplices) that successively pass through the n + 1 units placed on the n-simplex. Consequently, to model the n-order interactions arising from sequential participation, we need to analyze both the faces of the simplices (i.e. to ensure that all related units are indispensable) and the paths defined on these faces (i.e. to describe the sequential participation).

To describe the interactions propagating along the paths in n-simplices, we develop an operator named as the high-order path Laplacian $\mathbf{L}_{H}^{\left(n\right)}$

$$\begin{align} \left[\mathbf{L}_{H}^{\left(n\right)}\right]_{ij} = \frac{1}{n}\left(\delta_{ij}\mathbf{P}_{i}^{\left(n\right)}-\mathbf{B}_{ij}^{\left(n\right)}\right).\end{align} \tag{ 5 }$$

In equation (5), the coefficient $\frac{1}{n}$ denotes the normalization of the propagation speed of interactions along paths. The term $\mathbf{P}_{i}^{\left(n\right)}$ counts the number of paths that traverse all n + 1 units on the n-simplex without repetition and terminate at unit i

$$\begin{align} \mathbf{P}_{i}^{\left(n\right)} = n!\sum_{\left\{i\right\}_{n+1}\subset V}\prod_{\omega\in S^{\left\{i\right\}}_{n+1}}\mathbf{M}_{\omega},\end{align} \tag{ 6 }$$

where $n!$ counts the possibilities for the remaining part of the units in $\{i\}_{n+1}$ to form unique sequences after we fix unit v_i as the endpoint of these sequences. The term $\mathbf{B}_{ij}^{\left(n\right)}$ counts the number of paths that traverse all n + 1 units on the n-simplex without repetition, initiate at unit j, and terminate at unit i

$$\begin{align} \mathbf{B}_{ij}^{\left(n\right)} = \left(1-\delta_{ij}\right)\left(n-1\right)!\sum_{\left\{i,j\right\}_{n+1}\subset V}\prod_{\omega\in S^{\left\{i,j\right\}}_{n+1}}\mathbf{M}_{\omega},\end{align} \tag{ 7 }$$

where the coefficient $\left(n-1\right)!$ is derived from the fact that the remaining part of the units in $\{i,j\}_{n+1}$ can form $\left(n-1\right)!$ unique sequences if we fix units v_i and v_j as the endpoints.

The construction process of the high-order path Laplacian $\mathbf{L}_{H}^{\left(n\right)}$ is elaborated in section II in the supplementary materials. See figure 1(e) for several instances of the high-order path Laplacian. In a special case with n = 1, equations (6) and (7) reduce to the classic degree vector and adjacency matrix, respectively (i.e. we have $\mathbf{P}^{\left(1\right)} = \mathbf{D}^{\left(1\right)}$ and $\mathbf{B}^{\left(1\right)} = \mathbf{A}$). Therefore, equation (5) denotes the classic Laplacian when interactions are pairwise. In most general cases with n > 1, the n-order interactions formed by sequential participation are not equivalent to those caused by simultaneous participation (see section II in the supplementary materials for details).

In general, one can choose the multi-order Laplacian $\mathbf{L}_{M}^{\left(n\right)}$ and the high-order path Laplacian $\mathbf{L}_{H}^{\left(n\right)}$ according to application demands. As we have suggested above, while $\mathbf{L}_{M}^{\left(n\right)}$ conveys information about the high-order interactions that units are simultaneously engaged with, the definition of $\mathbf{L}_{H}^{\left(n\right)}$ is more applicable to the cases where a high-order interaction is formed by an irreducible sequence of the participation of related units. In section II in the supplementary materials, we explain how these two operators differ from each other in governing the dynamics of system evolution (e.g. with different decay rates of different orders). Specifically, the decay rate of the system evolution of the simultaneous participation (i.e. $\mathbf{L}_{M}^{\left(n\right)}$) is same as that of the sequential participation (i.e. $\mathbf{L}_{H}^{\left(n\right)}$) when interactions are pairwise (i.e. n = 1). When units exhibit two-order interactions, the sequential participation implies slower decay rates of the system evolution compared with the simultaneous participation. When interactions occur at higher orders (i.e. n > 2), the sequential participation leads to a larger decay rate. These differences highlight the dissimilar physics meanings and characteristics of $\mathbf{L}_{M}^{\left(n\right)}$ and $\mathbf{L}_{H}^{\left(n\right)}$.

Based on the unified mathematical forms presented in sections 2.1 and 2.2, we have developed efficient code implementations of these two kinds of operators, which can be seen in [51]. Meanwhile, in section III in the supplementary materials, we explain the similarities and differences between our considered Laplacian operators and other Laplacian operators used in topology theories (e.g. the combinatorial Laplacian [60] and the Hodge Laplacian [23, 61] in persistent homology analysis [62–67]), which helps us understand the theoretical significance of designing an RG on $\mathbf{L}_{M}^{\left(n\right)}$ and $\mathbf{L}_{H}^{\left(n\right)}$. For convenience, we uniformly denote them as $\mathbf{L}^{\left(n\right)}$ in our subsequent derivations. One can specify the actual definition of $\mathbf{L}^{\left(n\right)}$ in the application.

To offer a clear vision, we first introduce our approach, referred to as the SRG, in an ergodic case, as assumed by [10]. The SRG renormalizes a system of the p-order according to the properties of q-order interactions ($p\unicode{x2A7D} q$). Specifically, the renormalization is realized as the following:

• (1)  
  In each kth iteration, we have two iterative Laplacian operators, $\mathbf{L}^{\left(p\right)}_{k}$ and $\mathbf{L}^{\left(q\right)}_{k}$, and their associated high-order network sketches, $\mathbf{G}^{\left(p\right)}_{k}$ and $\mathbf{G}^{\left(q\right)}_{k}$, with $N^{\left(p\right)}_{k} = N^{\left(q\right)}_{k}$ units (there is an edge connecting between two units in $\mathbf{G}^{\left(p\right)}_{k}$ or $\mathbf{G}^{\left(q\right)}_{k}$ only if they share a p- or q-order interaction).
• (2)  
  In the real space, we first search for targets to coarse grain. Specifically, we search for every cluster of units sharing short-range high-order interactions within the timescale $\widehat{\tau}^{\left(q\right)}$, such that we can coarse grain these units into a macro-unit. Hence, we generate a reference network, $\mathbf{H}^{\left(q\right)}_{k}$, as follows. We first initialize $\mathbf{H}^{\left(q\right)}_{k}$ as a null network (i.e. with no edge). We follow equation (11) to calculate the density operator associated with the timescale $\widehat{\tau}^{\left(q\right)}$ as $\widehat{\rho}^{\,{\left(q\right)}}_{k} = \exp(-\widehat{\tau}^{\left(q\right)}\mathbf{L}^{\left(q\right)}_{k})/\operatorname{tr}(\exp(-\widehat{\tau}^{\left(q\right)}\mathbf{L}^{\left(q\right)}_{k}))$. In operator $\widehat{\rho}^{\,{\left(q\right)}}_{k}$, its main diagonal elements, related to self-interactions of units, are frequently greater than most of its off-diagonal elements (e.g. see figure 2(b)). If the $\left(i,j\right)$th element is greater than the $\left(i,i\right)$th or the $\left(j,j\right)$th element, then the q-order interactions between units v_i and v_j are sufficiently quick such that their accumulations within the timescale $\widehat{\tau}^{\left(q\right)}$ are faster than those of self-interactions. In this case, we add an edge between units v_i and v_j in network $\mathbf{H}^{\left(q\right)}_{k}$. By progressively implementing the edge-adding procedure for every pair of units, we can form $n_{k}^{\left(q\right)}$ clusters in $\mathbf{H}^{\left(q\right)}_{k}$. Each cluster contains a set of strongly correlated units with short-range high-order interactions. Basically, these clusters can be understood as Kadanoff blocks (i.e. the quasi-independent blocks, whose mean size serves as the correlation length of the system) in the real space renormalization procedure [71–73].
• (3)  
  After we find the targets for coarse graining in the real space, we implement a renormalization procedure. Specifically, in $\mathbf{G}^{\left(q\right)}_{k}$, each set of units that belong to the same cluster in $\mathbf{H}^{\left(q\right)}_{k}$ are aggregated into a macro-unit to generate a new network $\mathbf{G}^{\left(q\right)}_{k+1}$. Consequently, $n^{\left(q\right)}_{k}$ macro-units remain in $\mathbf{G}^{\left(q\right)}_{k+1}$ after coarse graining. In $\mathbf{G}^{\left(q\right)}_{k+1}$, two macro-units, v_i and v_j , are connected if at least one unit aggregated into v_i is connected with one unit aggregated into v_j in $\mathbf{G}^{\left(q\right)}_{k}$. If a unit or an interaction in $\mathbf{G}^{\left(q\right)}_{k}$ does not participate in any coarse graining according to $\mathbf{H}^{\left(q\right)}_{k}$, it is adopted into $\mathbf{G}^{\left(q\right)}_{k+1}$ without modification. Parallel to this process, we also follow the same rules to aggregate the units in $\mathbf{G}^{\left(p\right)}_{k}$ that belong to the same cluster of $\mathbf{H}^{\left(q\right)}_{k}$ into a macro-unit in $\mathbf{G}^{\left(p\right)}_{k+1}$. The generated $\mathbf{G}^{\left(p\right)}_{k+1}$ contains $n^{\left(q\right)}_{k}$ macro-units as well. Therefore, we define $N^{\left(p\right)}_{k+1} = N^{\left(q\right)}_{k+1} = n^{\left(q\right)}_{k}$.
• (4)  
  In the moment space, we first search for modes to reduce. Specifically, we find that the eigenvalues of $\mathbf{L}^{\left(q\right)}_{k}$ that are smaller than $\frac{1}{\widehat{\tau}^{\left(q\right)}}$, where $\widehat{\tau}^{\left(q\right)}$ is the timescale for q-order simplicial complexes (see section 4.2 for the selection of $\widehat{\tau}^{\left(q\right)}$). These eigenvalues correspond to long-range high-order interactions because the contributions of their eigenvectors have slow decay rates in equation (10). We denote the number of these eigenvalues as $s_{k}^{\left(q\right)}$. Other eigenvalues correspond to short-range high-order interactions and can be coarse grained. Moreover, according to equation (11), the totality of long-range high-order interactions (i.e. corresponding to smaller eigenvalues) is generally greater than that of short-range ones (i.e. associated with larger eigenvalues). Therefore, integrating out short-range high-order interactions does not damage the effective information contained in the system significantly. The only problem is that the eigenvalue spectrum cutoff defined above may overestimate the number of short-range high-order interactions and make $s_{k}^{\left(q\right)}$ too small. This is because the eigenvalue spectrum of $\mathbf{L}^{\left(q\right)}_{k}$ does not reflect the possibility that the short-range q-order interactions between units v_i and v_j can be slower than their self-interactions. For instance, although the interactions between units v_i and v_j are classified as short-range, the $\left(i,j\right)$th element of $\widehat{\rho}^{\,{\left(q\right)}}_{k}$ can be smaller than both the $\left(i,i\right)$th and the $\left(j,j\right)$th elements. In this case, the timescale $\widehat{\tau}^{\left(q\right)}$ may be too short for the correlations between v_i and v_j to accumulate to a non-negligible level. Compared with self-interactions, the interactions between units v_i and v_j are less important in determining the behaviors of v_i and v_j since these interactions are not strong enough. To avoid the possible overestimation, we use $n_{k}^{\left(q\right)}$, the number of Kadanoff blocks in $\mathbf{H}^{\left(q\right)}_{k}$, as a bound. In the case with overestimation, the number of long-range interactions is underestimated (i.e. $s_{k}^{\left(q\right)}\unicode{x2A7D} n_{k}^{\left(q\right)}$) and we correct it by defining $m_{k}^{\left(q\right)} = \max\{s_{k}^{\left(q\right)},n_{k}^{\left(q\right)}\}$. The value of $m_{k}^{\left(q\right)}$ defines the number of modes to keep during renormalization.
• (5)  
  After we find the modes to reduce in the moment space, we exclude these modes with fast decay rates from the system such that $m_{k}^{\left(q\right)}$ modes remain. Specifically, the current q-order Laplacian, $\mathbf{L}^{\left(q\right)}_{k}$, is reduced to the rescaled contributions of the eigenvectors associated with long-range high-order interactions

  $$\begin{align} \mathbf{Q}^{\left(q\right)}_{k} = \sum_{\lambda^{\left(q\right)}_{k}}I_{\Lambda^{\left(q\right)}_{k}}\left(\lambda^{\left(q\right)}_{k}\right)\lambda^{\left(q\right)}_{k}\big\vert\lambda^{\left(q\right)}_{k}\big\rangle\big\langle\lambda^{\left(q\right)}_{k}\big\vert,\end{align} \tag{ 17 }$$

  where $\Lambda^{\left(q\right)}_{k}$ is the set of the $m_{k}^{\left(q\right)}$ smallest eigenvalues of $\mathbf{L}^{\left(q\right)}_{k}$ (i.e. the first $m_{k}^{\left(q\right)}$ eigenvalues after we sort the eigenvalue spectrum in a decreasing order), and $I_{A}\left(\cdot\right)$ is the indicator function defined on an arbitrary set A. Meanwhile, the current p-order Laplacian, $\mathbf{L}^{\left(p\right)}_{k}$, is reduced to the contributions of the $m_{k}^{\left(q\right)}$ smallest eigenvalues (i.e. long-range p-order interactions)

  $$\begin{align} \mathbf{Q}^{\left(p\right)}_{k} = \sum_{\omega^{\left(p\right)}_{k}}I_{\Omega^{\left(p\right)}_{k}}\left(\omega^{\left(p\right)}_{k}\right)\omega^{\left(p\right)}_{k}\big\vert\omega^{\left(p\right)}_{k}\big\rangle\big\langle\omega^{\left(p\right)}_{k}\big\vert,\end{align} \tag{ 18 }$$

  where $\Omega^{\left(p\right)}_{k}$ is the set of the $m_{k}^{\left(q\right)}$ smallest eigenvalues of $\mathbf{L}^{\left(p\right)}_{k}$. Note that the eigenvalues in $\Omega^{\left(p\right)}_{k}$ are not necessarily smaller than $\frac{1}{\widehat{\tau}^{\left(q\right)}}$.
• (6)  
  The correspondence between moment space and real space for q-order interactions is ensured by a Laplacian $\mathbf{L}^{\left(q\right)}_{k+1}$. For v_i and v_j , two macro-units in $\mathbf{G}^{\left(q\right)}_{k+1}$, we define

  $$\begin{align} \left[\mathbf{L}^{\left(q\right)}_{k+1}\right]_{v_{i}v_{j}} = \big\langle v_{i}\big\vert\mathbf{Q}^{\left(q\right)}_{k}\big\vert v_{j}\big\rangle,\end{align} \tag{ 19 }$$

  if v_i and v_j are connected in $\mathbf{G}^{\left(q\right)}_{k+1}$. We set $[\mathbf{L}^{\left(q\right)}_{k+1}]_{v_{i}v_{j}} = 0$ if v_i and v_j are disconnected in $\mathbf{G}^{\left(q\right)}_{k+1}$. Here, $\big\vert v_{j}\big\rangle$ is a $N^{\left(q\right)}_{k}$-dimensional ket, where unitary components correspond to all the units in $\mathbf{G}^{\left(q\right)}_{k}$ that are aggregated into v_i and zero components correspond to all other units in $\mathbf{G}^{\left(q\right)}_{k+1}$ that are not covered by ν_i (i.e. the value of $[\mathbf{L}^{\left(q\right)}_{k+1}]_{v_{i}v_{j}}$ is the sum of all $[\mathbf{Q}^{\left(q\right)}_{k}]_{xy}$ if v_i and v_j are connected, where x goes through all the units in $\mathbf{G}^{\left(q\right)}_{k}$ that are aggregated into v_i and y goes through all the units in $\mathbf{G}^{\left(q\right)}_{k}$ that are aggregated into v_j ). After defining all off-diagonal elements, the diagonal elements of $\mathbf{L}^{\left(q\right)}_{k+1}$ are defined by following

  $$\begin{align} \left[\mathbf{L}^{\left(q\right)}_{k+1}\right]_{v_{i}v_{i}} = -\sum_{v_{j}}\left[\mathbf{L}^{\left(q\right)}_{k+1}\right]_{v_{i}v_{j}}.\end{align} \tag{ 20 }$$

  Mathematically, this procedure actually defines a similarity transformation T between $\mathbf{Q}^{\left(q\right)}_{k}$ and $\mathbf{L}^{\left(q\right)}_{k+1}$

  $$\begin{align} T^{-1}\mathbf{Q}^{\left(q\right)}_{k}T = \operatorname{diag}\left(\left[\mathbf{L}^{\left(q\right)}_{k+1},\mathbf{0}\right]\right),\end{align} \tag{ 21 }$$

  where 0 is a $(N^{\left(q\right)}_{k}-n^{\left(q\right)}_{k})$-dimensional all-zero square matrix. It is clear that the similarity transformation is given as $T = (\big\vert v_{1}\big\rangle,\ldots, \big\vert v_{n^{\left(q\right)}_{k}}\big\rangle, \big\vert \eta_{n^{\left(q\right)}_{k}+1}\big\rangle,\ldots, \big\vert \eta_{N^{\left(q\right)}_{k}}\big\rangle)$, where each ket η_j is selected such that the columns of T define a group of orthonormal bases. Similarly, we can derive the Laplacian $\mathbf{L}^{\left(p\right)}_{k+1}$ for p-order interactions based on $\mathbf{Q}^{\left(p\right)}_{k}$. The derived $\mathbf{L}^{\left(p\right)}_{k+1}$ and $\mathbf{L}^{\left(q\right)}_{k+1}$ are used in the $\left(k+1\right)$th iteration.

The renormalization of p-order interactions in both the real space (i.e. see steps (2) and (3)) and the moment space (i.e. see steps (4) and (5)) is realized under the guidance of q-order interactions. After renormalizing $\mathbf{L}^{\left(p\right)}_{k}$ by repeating steps (1)–(6), the SRG progressively drives the system toward an intrinsic scale of p-order interactions that exceeds the microscopic scales. More specifically, the SRG creates a q-order-guided iterative scale transformation that, irrespective of starting from which local scale, always leads to the latent critical point associated with p-order interactions (if there is any). It provides opportunities to verify whether the concerned thermodynamic functions are scale-invariant or not. See figure 3(a) for a summary of key steps in our approach. Note that the use of a reference network in step (2) is not necessary for theoretical derivations, yet it is favorable for computer programming.

Zoom In Zoom Out Reset image size

In a non-ergodic case, the presented procedure cannot be directly adopted. This is because our proposed renormalization approach is rooted in the relations among the time evolution operators defined by equations (8)–(10), the exponential decay of diffusion modes in the moment space, and the evolution paths in the path space, which implicitly depend on the assumption of the ergodicity of system states. In the non-ergodic case where a set of states can never evolve to or be transformed from another set of states (i.e. the network connectivity is absent), path integrals between these state sets are ill-posed and the diffusion over the network can be subdivided into several irrelevant sub-processes (i.e. the diffusion processes defined on two disconnected clusters are independent from each other). Dealing with these properties is crucial for our work since the SRG is oriented toward analyzing high-order interactions, which are usually sparsely distributed and may lack ergodicity.

Here, we develop a divide-and-conquer approach. Our idea arises from the fact that we can treat the diffusion from one cluster (i.e. a connected component) to another cluster on any order as being decelerated to a condition with zero velocity (i.e. being infinitely slow). This infinitely long-term diffusion process contains no short-range high-order interactions and, therefore, should not be integrated out during renormalization. Consequently, the SRG does not need to deal with the diffusion across any pair of clusters and can focus on the diffusion within each cluster severally. In the k-iteration, we consider the following two cases:

• (A)  
  If $\mathbf{G}^{\left(q\right)}_{k}$ is connected (i.e. the ergodicity holds), we directly apply steps (1)–(6) on $\mathbf{G}^{\left(q\right)}_{k}$ to guide the renormalization of p-order interactions and derive the inputs of the $\left(k+1\right)$-iteration.
• (B)  
  If $\mathbf{G}^{\left(q\right)}_{k}$ is disconnected (i.e. the ergodicity does not hold) and has r clusters denoted by $C_{1},\;\ldots,C_{r}$, we respectively apply steps (1)–(6) on each cluster of $\mathbf{G}^{\left(q\right)}_{k}$ to guide the renormalization of p-order interactions formed among the units that belong to this cluster. Specifically, we treat each cluster C_i as a network and input it into steps (1)-(6) with $\mathbf{G}^{\left(p\right)}_{k}$. Note that the above procedure does not affect $\mathbf{G}^{\left(p\right)}_{k}$ if the selected cluster C_i is trivial (i.e. contains only one unit). After dealing with all r clusters, the obtained results are used for the $\left(k+1\right)$-iteration.

In figure 3(b), we conceptually illustrate our approach, which enables us to relax the constraint on system ergodicity while implementing the SRG.

In section III in the supplementary materials, we explain how the SRG is related to the persistent homology theory [62, 74–76] and suggest its potential applicability to analyzing the lifetime of topological properties. In section IV in the supplementary materials, we discuss the SRG in the form of a conventional RG, suggesting the benefits of including high-order interactions into RG analysis, and indicate the relations between the SRG and other frameworks [31, 77–79]. Specifically, we find that the ideas underlying the SRG can be clearly interpreted by relating the SRG with a PCA-like RG theory [31].

After developing the renormalization procedure, we attempt to present a detailed analysis of the timescale $\tau^{\left(q\right)}$ and its associated scaling relation.

In the ergodic case, the setting of $\tau^{\left(q\right)}$ in the SRG (as we have mentioned in step (1)) is finished in the first iteration and required to make the specific heat, an indicator of the transitions of diffusion processes, maximize (see similar ideas in [10]). In section IV in the supplementary materials, we explain why it is valid to define renormalization flows based on the specific heat. Specifically, we consider the q-order spectral entropy in the first iteration of the SRG (see [42] for the definition of spectral entropy)

$$\begin{align} S^{\left(q\right)}_{1}\left(\tau^{\left(q\right)}\right) & =-\operatorname{tr}\left[\rho^{\left(q\right)}_{1}\log \left(\rho^{\left(q\right)}_{1}\right)\right],\end{align} \tag{ 22 }$$

$$\begin{align} & = \log \left(\sum_{\lambda^{\left(q\right)}_{1}}\exp\left(-\tau^{\left(q\right)}\lambda^{\left(q\right)}_{1}\right)\right)+\tau^{\left(q\right)}\left\langle\lambda^{\left(q\right)}_{1}\right\rangle_{\rho},\end{align} \tag{ 23 }$$

where each $\lambda^{\left(q\right)}_{1}$ is an eigenvalue of $\mathbf{L}^{\left(q\right)}_{1}$, and $\rho^{\left(q\right)}_{1}$ is the density operator derived from $\mathbf{L}^{\left(q\right)}_{1}$ following equation (11). Note that we have $\langle\lambda^{\left(q\right)}_{1}\rangle_{\rho} = \operatorname{tr}(\mathbf{L}^{\left(q\right)}_{1}\rho^{\left(q\right)}_{1})$. The q-order specific heat is defined according to the first derivative of the q-order spectral entropy

$$\begin{align} X^{\left(q\right)}_{1}\left(\tau^{\left(q\right)}\right)& =-\frac{\mathsf{d}}{\mathsf{d}\log\left(\tau^{\left(q\right)}\right)}S^{\left(q\right)}_{1}\left(\tau^{\left(q\right)}\right),\end{align} \tag{ 24 }$$

$$\begin{align} & = -\left(\tau^{\left(q\right)}\right)^{2}\frac{\mathsf{d}}{\mathsf{d}\tau^{\left(q\right)}}\left\langle\lambda^{\left(q\right)}_{1}\right\rangle_{\rho}.\end{align} \tag{ 25 }$$

See section V in the supplementary materials for detailed derivations. The value of $\tau^{\left(q\right)}$ is selected to ensure $\frac{\mathsf{d}}{\mathsf{d}\tau^{\left(q\right)}}X^{\left(q\right)}_{1}\left(\tau^{\left(q\right)}\right) = 0$, which indicates an infinite deceleration condition of diffusion processes. See figure 3(c) for illustrations. Apart from the infinite deceleration condition, determination of the timescale value $\widehat{\tau}^{\left(q\right)}$ used in step (2) of the SRG also requires us to consider the high-order properties of diffusion. For convenience, we denote $\widetilde{\tau}^{\left(q\right)}$ as the value of $\tau^{\left(q\right)}$ under the infinite deceleration condition (i.e. there is $\frac{\mathsf{d}}{\mathsf{d}\tau^{\left(q\right)}}X^{\left(q\right)}_{1}\left(\tau^{\left(q\right)}\right)\big\vert_{\widetilde{\tau}^{\left(q\right)}} = 0$). The optimal timescale for the SRG is set as

$$\begin{align} \widehat{\tau}^{\left(q\right)} = \frac{2\varepsilon}{q+1}\widetilde{\tau}^{\left(q\right)},\end{align} \tag{ 26 }$$

which is the average subdivided timescale for each 1-simplex in an n-simplex (i.e. the average number of adjacent 1-simplices between arbitrary units i and j is $\frac{n+1}{2}$). Note that $\widehat{\tau}^{\left(q\right)} = \widetilde{\tau}^{\left(q\right)}$ holds for pairwise interactions. The value of ε is used to correct the errors of timescale calculation when system units share ultra-dense interactions, which arise from the numerical bias of diffusion rate estimation in ultra-dense systems. The definition of ε can be seen in section VI in the supplementary materials.

The above procedure can be directly generalized to the non-ergodic case, where we treat each cluster C_i as a network and derive a $\widehat{\tau}^{\left(q\right)}$ for it if it has not been assigned with a timescale yet. The timescale of a cluster does not change after it is assigned. See figure 3(c) for illustrations.

After determining the timescale, the SRG proposed in section 4.1 can be applied to renormalize interacting systems in a non-parametric manner. See section XIII in the supplementary materials for the programmatic implementation of the SRG.

In figure 4(a), we explain how the SRG iterates. For an input system, we define p and q to realize the renormalization guided by different orders of interactions. Given each pair of $\left(p,q\right)$, we first search p- and q-simplices to represent the system of the p- and q-order (i.e. generate the high-order network sketches, $\mathbf{G}^{\left(p\right)}_{1}$ and $\mathbf{G}^{\left(q\right)}_{1}$, as shown in figure 1). Meanwhile, we derive the critical timescale, $\widehat{\tau}^{\left(q\right)}$, on $\mathbf{G}^{\left(q\right)}_{1}$ following section 4.2. After the critical timescale is calculated, we iteratively realize steps (1)–(6) introduced in section 4.1 to renormalize the system once, twice, and so on. For instance, in the first iteration of the SRG, we follow steps (2)–(3) to find the Kadanoff blocks of the q-order (i.e. clusters in $\mathbf{H}^{\left(q\right)}_{1}$) and coarse grain the units in each Kadanoff block into a macro-unit on the p and q-order to derive $\mathbf{G}^{\left(p\right)}_{2}$ and $\mathbf{G}^{\left(q\right)}_{2}$. Then, we follow steps (4)–(6) to calculate the associated high-order Laplacian operators $\mathbf{L}^{\left(p\right)}_{2}$ and $\mathbf{L}^{\left(q\right)}_{2}$. The obtained $\mathbf{G}^{\left(p\right)}_{2}$, $\mathbf{G}^{\left(q\right)}_{2}$, $\mathbf{L}^{\left(p\right)}_{2}$, and $\mathbf{L}^{\left(q\right)}_{2}$ serve as inputs for the SRG to run the second iteration to generate $\mathbf{G}^{\left(p\right)}_{3}$, $\mathbf{G}^{\left(q\right)}_{3}$, $\mathbf{L}^{\left(p\right)}_{3}$, and $\mathbf{L}^{\left(q\right)}_{3}$. In this manner, the iteration continues and a renormalization flow is gradually generated.

Zoom In Zoom Out Reset image size

In figures 4(b)–(d), we mainly illustrate how the real space coarse graining in the SRG occurs at different orders. We define p = 1 and $q\in\{1,2,3\}$ to generate renormalization flows and visualize the changes caused by coarse graining the units within each Kadanoff block. Two important messages are conveyed by figure 4. First, whenever a Kadanoff block of m units is found in $\mathbf{H}^{\left(q\right)}_{k}$, the size of $\mathbf{G}^{\left(p\right)}_{k+1}$ is reduced by m − 1 (this property also holds for $\mathbf{G}^{\left(q\right)}_{k+1}$). If there is no Kadanoff block in $\mathbf{H}^{\left(q\right)}_{k}$ (i.e. no unit to coarse grain), the SRG arrives at its fixed point and the system remains unchanged during subsequent renormalization. Second, as shown in the comparison across different values of q, the system reduction during renormalization is always smaller given a larger value of q, which mainly arises from the sparser distribution of high-order interactions.

To demonstrate the meaning of distinguishing between p and q during renormalization, we compare between the cases with p ≠ q (i.e. renormalize the system of an order according to its properties in another order) and p = q (i.e. equivalent to the situation where we do not distinguish between p and q) in section VII in the supplementary materials. In Sfigure 2 in the supplementary materials, we observe significant differences in both real space (i.e. $\mathbf{G}^{\left(p\right)}_{k}$) and moment space (i.e. $\mathbf{L}^{\left(p\right)}_{k}$) between p = q and p ≠ q, demonstrating the effectiveness of supporting p ≠ q during renormalization. Therefore, by distinguishing between p and q, the SRG is enabled to analyze more phenomena than classic frameworks with $p\equiv q$. Meanwhile, these results qualitatively suggest the effects of higher-order interactions on lower-order interactions during renormalization (i.e. in the case with p ≠ q), whose quantification is explored in our subsequent analysis.

In figure 5, we present more instances of the renormalization flow generated by the SRG. We apply the SRG to synthetic interacting systems whose pairwise interactions follow the random tree (RT), the square lattice (SL), the triangular lattice (TL), the Watts–Strogatz networks (WS), and the Erdos–Renyi networks (ER). Among these random networks, the SL and TL can be defined with or without the periodic boundary conditions (PBCs). See section I in the supplementary materials for the definitions of these random networks. As shown in figure 5, the key structural properties of the RT, SL, and TL are principally preserved during renormalization (i.e. we can see qualitative similarities in $\mathbf{G}^{\left(p\right)}_{k}$ across different iterations) while the SRG implemented on the WS and ER networks gradually reveals specific latent structures that are dissimilar to the original ones. Among these instances, the lattice systems with PBCs are observed to make the renormalization ineffective. This is a phenomenon caused by the identical properties of all units under the PBCs (i.e. all units are exactly the same as each other at any q-order). Because the renormalization procedure in the SRG begins by distinguishing between short- and long-range interactions, the SRG can effectively reduce a system only when units and interactions are not identical. Otherwise, it is impossible to tell which interactions are long-range and which are short-range. In figures 5(c) and (e), the SRG does not coarse grain systems because it discovers no short-range interaction. This phenomenon does not cause conflicts in our work because a lattice system with identical units and interactions can be treated as an extreme case of self-similarity to a certain extent (i.e. each sub-system is identical to the whole system) and can be invariant for coarse graining. To avoid this trivial phenomenon in our analysis, all the lattice systems implemented in our subsequent experiments (e.g. figures 6 and 7) are defined without the PBCs.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Given the SRG, a crucial task is to analyze the potential scaling relation and quantify the effects of high-order interactions on it.

In the ergodic case, we denote $\widetilde{X}^{\left(q\right)}$ and $\widetilde{\rho}^{\left(q\right)}$ as the specific heat and density operator associated with $\widetilde{\tau}^{\left(q\right)}$. We can insert these variables into equation (25) to derive the following relation

$$\begin{align} \frac{\widetilde{X}^{\left(q\right)}}{\widetilde{\tau}^{\left(q\right)}} = \left\langle\lambda^{\left(q\right)}_{1}\right\rangle_{\widetilde{\rho}}+\mu,\;\forall \mu\in\mathbb{R},\end{align} \tag{ 27 }$$

where $\langle\lambda^{\left(q\right)}_{1}\rangle_{\widetilde{\rho}} = \operatorname{tr}(\mathbf{L}^{\left(q\right)}_{1}\widetilde{\rho}^{\left(q\right)})$ and the constant µ is obtained when we solve equation (25) as a differential equation. If the Laplacian eigenvalue spectrum (i.e. the density of the eigenvalues of $\mathbf{L}^{\left(q\right)}_{1}$) follows a general power-law form, $\operatorname{Prob}(\lambda^{\left(q\right)}_{1})\sim (\lambda^{\left(q\right)}_{1})^{\beta^{\left(q\right)}}$, we can obtain a scaling relation

$$\begin{align} \beta^{\left(q\right)} = \widetilde{X}^{\left(q\right)}-1-\nu,\;\exists \nu\in\left[0,\infty\right),\end{align} \tag{ 28 }$$

whose derivations are presented in section VIII in the supplementary materials. The parameter ν is a specific constant. For a system with scale-invariance properties of the q-order, the power-law form of the Laplacian eigenvalue spectrum holds and should be invariant under the transformation of the SRG. Therefore, we expect to see an approximately constant $\widetilde{X}^{\left(q\right)}$ in a certain timescale interval (i.e. the infinite deceleration condition), accompanied by a fixed $\beta^{\left(q\right)}$ during renormalization. For a system without scale invariance, the scaling relation does not hold true.

To generalize our analysis to the non-ergodic case, we respectively derive $\widetilde{\tau}^{\left(q\right)}\left(i\right)$ and $\beta^{\left(q\right)}\left(i\right)$ for each cluster C_i of $\mathbf{G}^{\left(q\right)}_{1}$ by following equations (27) and (28). These definitions enable us to obtain

$$\begin{align} \left\langle\lambda^{\left(q\right)}_{1}\right\rangle_{\widetilde{\rho}} = \left\langle\frac{\beta^{\left(q\right)}\left(i\right)+1}{\widetilde{\tau}^{\left(q\right)}\left(i\right)}\right\rangle_{i}+\nu,\;\exists \nu\in\left[0,\infty\right),\end{align} \tag{ 29 }$$

and

$$\begin{align} \widetilde{X}^{\left(q\right)} = \widetilde{\tau}^{\left(q\right)}\left\langle\frac{\beta^{\left(q\right)}\left(i\right)+1}{\widetilde{\tau}^{\left(q\right)}\left(i\right)}\right\rangle_{i}+\nu,\;\exists \nu\in\left[0,\infty\right),\end{align} \tag{ 30 }$$

where $\langle\cdot\rangle_{i}$ denotes weighted averaging across all clusters of $\mathbf{G}^{\left(q\right)}_{1}$ (i.e. weighted according to cluster size). In the non-ergodic case, variables $\widetilde{X}^{\left(q\right)}$ and $\widetilde{\tau}^{\left(q\right)}$ refer to the global specific heat and the timescale measured on $\mathbf{L}^{\left(q\right)}_{1}$ under the infinite deceleration conditions. Equations (27) and (28) actually serve as special cases of equations (29) and (30), where only one cluster exists. See section IX in the supplementary materials for derivations of these equations. Moreover, although the scaling relations in equations (27)–(30) are derived in a case with k = 1, we can also explore these scaling relations in $\mathbf{G}^{\left(q\right)}_{k}$ for k > 1 (i.e. we directly replace $\mathbf{G}^{\left(q\right)}_{1}$ with $\mathbf{G}^{\left(q\right)}_{k}$ during the calculation of equations (27)–(30)). In summary, we can verify if a system obeys scaling relations and further explore whether the SRG preserves scaling relations.

Equations (27)–(30) are derived for q-order interactions during renormalization and may divide from the original scaling relation of p-order interactions. These p-order interactions, initially corresponding to $\widetilde{X}^{\left(p\right)}$ and $\beta^{\left(p\right)}$, are renormalized according to $\widetilde{X}^{\left(q\right)}$ and $\beta^{\left(q\right)}$ in the SRG. In the case with p < q, we can measure the effects of q-order interactions on the scaling relation of p-order interactions as

$$\begin{align} \zeta\left(p,q\right) = \vert\widetilde{X}^{\left(p\right)}-\widetilde{X}^{\left(q\right)}\vert,\end{align} \tag{ 31 }$$

which is applicable to both ergodic and non-ergodic cases. Meanwhile, similar to the scaling relation analysis, the measurement of high-order effects can be implemented in any iteration of the SRG.

Overall, the proposed SRG offers an opportunity to verify the potential existence of scale invariance at different orders and classify interacting systems according to their high-order scaling properties.

We aim to verify the multi-order scale invariance (i.e. being invariant under the SRG transformation at different orders). Specifically, being scale-invariant at an order requires the system to be located at a certain fixed point of the renormalization flow (i.e. the system behaves similarly across all scales), leading to scale-independent characteristics of the system described by power-law scaling. Specifically, if we describe system states via a set of observables, their measurements are expected to maintain invariance across all the iterations of the SRG. Due to the internal complexity of the system's different orders, we may see a rich variety of fixed points in empirical data. In an opposite where the system lacks a fixed point, its behavior significantly changes across scales.

In section X in the supplementary materials, we implement our verification according to the behaviors of the Laplacian eigenvalue spectrum (the behavior of degree distributions is also presented as auxiliary information). Based on the existence of invariant power-law forms of Laplacian eigenvalue spectra under the SRG transformation, we can classify systems into scale-invariant, weakly scale-invariant, and scale-dependent types of different orders.

Apart from the above approach, there is a more direct way to realize the verification. Because a scale-invariant system of the q-order should invariably follow the scaling relation, being bound to satisfy equation (29) is one of the necessary conditions for being scale-invariant. Consequently, we can first measure the deviations from the scaling relation of the concerned systems before renormalization. Given significantly large deviations, the system cannot be scale-invariant, regardless of how it behaves during renormalization. In figure 6, the measurement is conducted on multiple types of systems in each kth iteration ($k\in\{1,2,3\}$) of the SRG (see section I in the supplementary materials for random network definitions, and see section XI in the supplementary materials for experimental details). As shown in figures 6(a)–(e), although the mean observations of $(\widetilde{X}^{\left(q\right)} ,\widetilde{\tau}^{\left(q\right)}\langle\frac{\beta^{\left(q\right)}\left(i\right)+1}{\widetilde{\tau}^{\left(q\right)}\left(i\right)}\rangle_{i}+\nu)$ averaged across all replicas collapse onto the scaling relation at every order and in each kth iteration, the departures from the scaling relation measured by the standard deviations of $\widetilde{X}^{\left(q\right)}$ and $\widetilde{\tau}^{\left(q\right)}\langle\frac{\beta^{\left(q\right)}\left(i\right)+1}{\widetilde{\tau}^{\left(q\right)}\left(i\right)}\rangle_{i}+\nu$ in some systems can be non-negligible. Specifically, as suggested in figure 6(f), the systems whose pairwise interactions follow known scale-invariant structures, such as the TL or the RT, do ensure that all their replicas follow the scaling relation with vanishing deviations in each kth iteration. The replicas of the systems whose pairwise interactions follow the Barabási–Albert network, a structure with weak scale-invariant properties (i.e. the Barabási–Albert network can be renormalized within a specific set of scales, but is not strictly scale-invariant [10]), have relatively small deviations at the first three orders as well. The replicas of other kinds of systems whose pairwise interactions follow the WS or the ER networks [80] exhibit larger deviations from the scaling relation at most orders, suggesting the absence of multi-order scale invariance.

In summary, the method reported in figure 6 can offer convenient detection of scale-dependent systems in practice because scale-invariant systems are expected to satisfy equation (29) before and during renormalization. To seek a more precise verification, one should apply the approach in section X in the supplementary materials.

Finally, we measure the high-order effects on the scaling relation by applying equation (31). As shown in figures 6(g)–(i), these high-order effects are non-negligible and generally increase with the difference between p and q. These results explain how the SRG differs from conventional RGs that focus on a single order.

In section III in the supplementary materials, we have explained how the SRG is theoretically related to the persistent homology theory and other concepts, where we suggest the potential applicability of the SRG to analyzing persistent homology. Below, we validate this applicability by using the SRG to explore the topological invariance.

In brief, the persistent homology theory can be used to study the lifetimes of different topological properties in a filtration process, during which we progressively add or remove simplices from the system [62, 74–76] (i.e. this is similar to the network evolution process in complex network studies [81]). Although the SRG is not exactly the same as the classic filtration approach because it modifies the system in a more complicated manner, it does offer an opportunity to study the evolution of topological properties during renormalization (i.e. we can measure different topological properties in $\mathbf{G}^{\left(p\right)}_{k}$ for each $k\unicode{x2A7E} 1$). Similar to the analysis of scale invariance, a topological property can be treated as invariant if it is persistent across all the iterations of the SRG (i.e. referred to as topological invariance). In our analysis, we mainly focus on the Betti numbers of different orders, a family of important topological properties used in biology [82–84] and physics [85–87] studies. In general, an m-order Betti number, $\mathcal{B}_{m}$, counts the number of m-dimensional holes in the system, which is a key topological characteristic of high-order structures and defines the dimension of mth homology [74]. In the terminology of graph theory, a 0-order Betti number counts the number of connected components and a 1-order Betti number is the number of cycles. Computationally, the m-order Betti number can be calculated using the toolbox offered by [88].

In figure 7, we explore the potential existence of topological invariance in different random network models, including the Barabási–Albert network (c = 4), the WS network (each unit initially has 10 neighbors and edges are rewired according to a probability of 0.3), the ER network (two units share an edge with a probability of 0.1), the TL, and the RT. See section I in the supplementary materials for their definitions. We generate 100 replicas for each random network model, where every replica consists of 200 units. In the experiment, we design an SRG with a certain combination of $\left(p,q\right)$ and let it function on each random network to generate a sequence of $\mathbf{G}^{\left(p\right)}_{k}$ for $k\in\{1,\ldots,5\}$. Given each $\mathbf{G}^{\left(p\right)}_{k}$, we measure the m-order Betti number, $\mathcal{B}_{m}$, on it ($m\in\{0,\ldots,3\}$). If an m-order Betti number is generally constant across all iterations, its corresponding topological properties can be treated as invariant during renormalization.

As shown in figure 7, the Betti numbers of scale-dependent systems (e.g. the WS and the ER networks) exhibit non-negligible changes on most orders when the SRG transforms systems across scales. Compared with these systems, weakly scale-invariant systems (e.g. the Barabási–Albert network) have less variable Betti numbers during renormalization. When systems become scale-invariant (e.g. the TL and the RT), the Betti numbers at all orders are generally invariant, except that some fluctuations may occur in certain replicas due to random effects. This finding is interesting because, to the best of our knowledge, there is no theoretical guarantee that the multi-order scale invariance is correlated with the topological invariance yet. According to our observations, we have the following speculation: specific topological invariance properties may exist near the critical point that makes the SRG, a framework dealing with high-order topological structures, have fewer effects on the information contained in the system. As a result, the behavior of the system seems to be approximately invariant under the scale transformation of the SRG, leading to the inference that the system is born at a fixed point of the renormalization flow. Exploring this speculation may deepen our understanding of the relation between the persistent homology [62, 74–76] and the criticality studied by the RG theory [6, 7], which remains a valuable direction for future works.

Next, we use the SRG to explore the roles of different orders of interactions in preserving or defining the organizational structures (e.g. latent communities) of complex systems, where we also validate the applicability of the SRG in organizational structure identification.

When organizational structures are given, we explore the conditions under which these structures can be maintained by the SRG. In section XII in the supplementary materials, we implement the verification using protein–protein interactions (e.g. orthologous and genetic interactions) in Caenorhabditis elegans [89]. We first extract the initial community structures of protein–protein interactions by applying the Louvain community detection algorithm [90] (see section I in the supplementary materials for an introduction). Then, we use the SRG to renormalize protein–protein interactions on a high order (i.e. with $\left(p,q\right) = \left(1,2\right)$) or a low order (i.e. with $\left(p,q\right) = \left(1,1\right)$). We define latent community structures according to unit aggregation during renormalization (see section XII in the supplementary materials for details). Consistent with [38], our results in section XII in the supplementary materials highlight that renormalization can never be treated as a trivial counterpart of clustering. Meanwhile, the SRG guided by high-order interactions (e.g. q = 2) can generally preserve the properties of initial community structures while the SRG guided by pairwise interactions does not. Based on these results, we speculate that high-order interactions are essential in forming and characterizing organizational structures. Communities cannot be treated as trivial collections of pairwise interactions, and it may be inappropriate to partition a system into sub-systems based only on pairwise relations. Given this observation, we naturally wonder if high-order interactions alone are sufficient to define ideal community structures (i.e. do not consider pairwise relations while determining system partition).

To tackle this question, we explore whether an SRG defined by different combinations of $\left(p,q\right)$ can be applied to identify latent communities when system partition is unknown. In figures 8(a)–(c), we define communities according to the renormalization flows of the SRG and compare them with the Louvain communities (i.e. the communities detected by the Louvain algorithm [90], see section I in the supplementary materials). Specifically, we classify the units into the same SRG community if they are aggregated into the same macro-unit during renormalization (e.g. we generate SRG communities after the first iteration of renormalization flows in figures 8(a)–(c)).

Zoom In Zoom Out Reset image size

To quantitatively compare the SRG communities with the Louvain communities, we analyze the correctness of communities (i.e. or referred to as the performance) [91] and the disconnectivity between communities (i.e. or referred to as the coverage) [91]. Note that both these concepts reflect the properties of communities on 1-order because they are defined on pairwise interactions [91] (see section I in the supplementary materials for precise definitions). As shown in figures 8(d)–(f), the behavior of the SRG is different from the Louvain community detection algorithm, which is designed to maximize the modularity of communities [90, 92, 93]. When the SRG is guided by pairwise interactions, the generated SRG communities tend to be disconnected from each other (i.e. lack pairwise interactions across communities) and have low accuracy (i.e. units in the same community may lack pairwise interactions because many short-range pairwise interactions are reduced). When the SRG is guided by high-order interactions, the generated SRG communities tend to be precise (i.e. units in the same community are always involved in high-order interactions) and non-isolated (i.e. across-community pairwise interactions exist since the SRG mainly reduces high-order interactions). We are inspired to generate the SRG communities by combining the properties of these two cases. We suggest considering pairwise interactions to correct the SRG communities generated by the SRG guided by high-order interactions. Specifically, for two SRG communities generated with $\left(p,q\right) = \left(1,2\right)$, if at least half of the units in one community have pairwise interactions with the units in another community, then all units in these two SRG communities are merged into the same community. The corrected SRG communities are presented in figures 8(d)–(f), exhibiting better capability in balancing between the correctness and the disconnectivity.

From another perspective, we quantify the consistency between the SRG communities and the Louvain communities in defining system partition in figures 8(g)–(i). As measured by the adjusted mutual information, completeness, and homogeneity (see section I in the supplementary materials for the definitions of these metrics, whose larger values suggest higher consistency extents), the partitions defined by the corrected SRG communities are generally consistent with, but not exactly equivalent to, those formed by the Louvain communities in all data sets. On the other hand, the SRG communities defined based only on high-order (i.e. with $\left(p,q\right) = \left(1,2\right)$) or low-order (i.e. with $\left(p,q\right) = \left(1,1\right)$) interactions do not necessarily maintain consistency with the Louvain communities in these data sets. Because the Louvain communities are defined by maximizing the modularity [90, 92, 93], we suggest that the optimization of modular structures cannot be realized based only on pairwise or high-order interactions. Instead, modular structures depend on multiple orders of interactions.

Taken together, high-order interactions are more effective than pairwise interactions in preserving the principal characteristics of organizational structures. However, high-order interactions alone are not sufficient to define ideal system partitions with optimal modular structures and clear separability among communities. For real systems, their organizational structures are defined by the intricate coexistence of various orders of interactions, exhibiting diverse behaviors across different orders. Although the SRG does not behave in a manner equivalent to the optimized algorithms for community detection (e.g. the Louvain algorithm [90]), it offers a flexible way to analyze the roles of the interactions of every order in characterizing organizational structures since we can freely select different combinations of $\left(p,q\right)$ in analysis. Moreover, after considering the interactions of multiple orders, the corrected SRG communities have the potential to offer a system partition scheme that is competitive with optimized algorithms (see figures 8(d)–(i)), which suggests the applicability of the SRG to identifying latent community structures of real systems.

Finally, we analyze the SRG in terms of informational properties. In previous studies, the analysis using information theory has suggested that RG frameworks essentially maximize the mutual information between relevant features and the environment [94, 95] or reduce the mutual information among irrelevant features [96]. In this work, we suggest exploring the SRG from the perspective of an information bottleneck [97] because an RG, similar to dimensionality reduction or representation learning approaches in machine learning [98–100], essentially deals with information encoding during data compression [101].

To realize this analysis, we follow one of our earlier works [102] to represent the associated high-order network sketch in the kth iteration of the renormalization flow, $\mathbf{G}^{\left(p\right)}_{k}$, by a Gaussian variable

$$\begin{align} \mathbf{X}^{\left(p\right)}_{k}\sim\mathcal{N}\left(\mathbf{0},\mathbf{L}^{\left(p\right)}_{k}+\frac{1}{N^{\left(p\right)}_{k}}\mathbf{J}\right),\end{align} \tag{ 32 }$$

where J is an all-one matrix and $N^{\left(p\right)}_{k}$ measures the number of units in $\mathbf{G}^{\left(p\right)}_{k}$. The derived Gaussian variable offers an optimal representation of $\mathbf{G}^{\left(p\right)}_{k}$ with network-topology-dependent smoothness and maximum entropy properties, and enables us to compare $\mathbf{G}^{\left(p\right)}_{k}$ across different iterations (see [102] for explanations). If necessary, one can further add a scaled unitary matrix, $a\mathbf{I}$ (a > 0), to the covariance matrix of $\mathbf{X}^{\left(p\right)}_{k}$ to ensure its semi-positive properties in the non-ergodic case (i.e. where $\mathbf{G}^{\left(p\right)}_{k}$ is disconnected).

We suggest that one considers the following information bottleneck

$$\begin{align} \mathcal{L}_{\textrm{IB}} = \underbrace{\mathcal{I}\left(\mathbf{X}^{\left(p\right)}_{k};\mathbf{X}^{\left(p\right)}_{1}\right)}_{\text{Preserved information}}-\psi\underbrace{\mathcal{I}\left(\mathbf{X}^{\left(p\right)}_{k};\mathbf{X}^{\left(p\right)}_{i}\right)}_{\text{Complexity}},\;\forall\psi \gt 0\end{align} \tag{ 33 }$$

where $i = k-1$ if k > 1 and i = k if k = 1. The parameter ψ defines the regularization strength. We denote $\mathcal{I}\left(\cdot;\cdot\right)$ as the mutual information (see section I in the supplementary materials for its definition), which can be estimated by the mutual information estimator with local non-uniformity corrections (we set the correction strength as 10⁻³ such that corrections only occur when necessary) [103]. This approach has been included in the non-parametric entropy estimation toolbox [104]. Note that equation (33) is different from the information bottleneck considered in [101] due to the discrepancy between our concerned question and [101].

Intuitively, the objective function in equation (33) evaluates whether the renormalized system in the kth iteration, $\mathbf{G}^{\left(p\right)}_{k}$, preserves the information of the original system, $\mathbf{G}^{\left(p\right)}_{1}$, by maximizing $\mathcal{I}(\mathbf{X}^{\left(p\right)}_{k};\mathbf{X}^{\left(p\right)}_{1})$ and reduces the trivial dependence on its previous state, $\mathbf{G}^{\left(p\right)}_{i}$, by minimizing $\mathcal{I}(\mathbf{X}^{\left(p\right)}_{k};\mathbf{X}^{\left(p\right)}_{i})$. For a valid RG, being able to preserve some information of the original system is an essential capability. Because a trivial solution of maximizing $\mathcal{I}(\mathbf{X}^{\left(p\right)}_{k};\mathbf{X}^{\left(p\right)}_{1})$ is to make $\mathbf{G}^{\left(p\right)}_{k} = \mathbf{G}^{\left(p\right)}_{k-1}$ for any k > 1 (i.e. there is no renormalization at all), we need to include the complexity term in equation (33) to avoid this trivial result.

Meanwhile, the possibility of maximizing $\mathcal{I}(\mathbf{X}^{\left(p\right)}_{k};\mathbf{X}^{\left(p\right)}_{1})$ depends on the properties of the system as well. If the system is scale-invariant, the preserved information does not rapidly vanish when we reduce the complexity since the information of $\mathbf{G}^{\left(p\right)}_{1}$ holds across scales (note that $\mathcal{I}(\mathbf{X}^{\left(p\right)}_{k};\mathbf{X}^{\left(p\right)}_{1})$ still reduces slightly because the decreasing dimensionality of $\mathbf{X}^{\left(p\right)}_{k}$ affects the numerical behavior of the estimator). If the system is scale-dependent, the preserved information crucially relies on the complexity and vanishes unless the complexity is high. Therefore, we can divide $\mathcal{I}(\mathbf{X}^{\left(p\right)}_{k};\mathbf{X}^{\left(p\right)}_{1})$ by $\mathcal{I}(\mathbf{X}^{\left(p\right)}_{k};\mathbf{X}^{\left(p\right)}_{i})$ to measure the complexity required by a single bit of preserved information in the kth iteration. A high value of complexity per information bit suggests that the system is scale-dependent because we cannot find a low-complexity coarse-grained system to represent it efficiently.

In figure 9(a), we show the information curves of the renormalization flows generated in four kinds of interacting systems. As suggested by our results, the renormalization flows of weakly scale-invariant (e.g. Barabási–Albert) and scale-invariant (e.g. RT) systems lead to sufficiently small complexity values and numerous preserved information bits, while the renormalization flows of scale-dependent systems (e.g. ER and WS) do not always achieve this (although the WS network realizes complexity reductions during renormalization, it is still not competitive with the Barabási–Albert network and RTs). The information bits of scale-dependent systems are preserved at the cost of maintaining high complexity, while scale-invariant and weakly scale-invariant systems progressively become more complexity-saving in preserving information during renormalization. These results can be quantitatively validated based on the complexity per information bit measured in figure 9(b). These findings suggest the possibility of verifying scale invariance from an informational perspective. Moreover, we observe that the SRGs guided by pairwise interactions (i.e. q = 1) frequently achieve more significant complexity reductions while preserving information. This phenomenon may arise from the fact that high-order interactions are distributed in systems in a sparser manner.

Zoom In Zoom Out Reset image size