On Complexity of Stability Analysis in Higher-order Ecological Networks through Tensor Decompositions

Abstract

Complex ecological networks are often characterized by intricate interactions that extend beyond pairwise relationships. Understanding the stability of higher-order ecological networks is salient for species coexistence, biodiversity, and community persistence. In this article, we present complexity analyses for determining the linear stability of higher-order ecological networks through tensor decompositions. We are interested in the higher-order generalized Lotka-Volterra model, which captures high-order interactions using tensors of varying orders. To efficiently compute Jacobian matrices and thus determine stability in large ecological networks, we exploit various tensor decompositions, including higher-order singular value decomposition, Canonical Polyadic decomposition, and tensor train decomposition, accompanied by in-depth computational and memory complexity analyses. We demonstrate the effectiveness of our framework with numerical examples.

I Introduction

Modeling complex systems and predicting their long-term behavior have drawn significant attention across diverse fields, including mathematics, physics, biology, and social science [29, 18]. Network representations, which assign individual components to nodes and their connections to edges, are widely used to characterize complex systems. They place a specific emphasis on understanding the influence of interactions, quantified by the edge weights, as in weighted Laplacian or adjacency matrices. Networks have been proven to be powerful tools in system modeling, with numerous works on related system-theoretic properties such as stability, controllability, and observability.

However, a comprehensive understanding of system characterizations with high-order interactions remains elusive. High-order interactions occur among a group of individual components, and thus can have a significant aggregated effect on the system, whereas standard networks only consider pairwise connections. To address this challenge, high-order networks (i.e., hypergraphs, see Fig. 1) and their tensorial representations are increasingly being used to model complex real-world systems, where adjacency coefficients can be defined on sets of nodes. Tensors, which are multidimensional arrays [22, 15, 26], have found applications across a wide array of fields, including dynamical systems [11, 13, 9], signal processing [27, 14], network science [10, 12, 31], data analysis [32], and more [30].

Not only can the static status of higher-order networks be studied, but also their long-standing interests in stability and state evolution. The associated dynamical system is thus composed of pairwise interaction matrix and high-order interaction tensors. Specifically, we are dealing with complex ecological networks, where the most commonly used model is known as the generalized Lotka–Volterra (GLV) model [7, 8]. It is defined as

$$\dot{\textbf{x}}=\textbf{x}*\Big{[}\textbf{r}+\textbf{A}\textbf{x}\Big{]},$$ (1)

where $\textbf{x}\in\mathbb{R}^{n}$ represents species abundance, $\textbf{r}\in\mathbb{R}^{n}$ is the vector of intrinsic growth rate, and $\textbf{A}\in\mathbb{R}^{n\times n}$ is the interaction matrix whose off-diagonal element $\textbf{A}_{ij}$ represents the effect of species $j$ upon species $i$. Throughout, we use $*$ to denote element-wise multiplication.

Numerous linear stability results have been developed for the GLV model by characterizing its community matrix (Jacobian matrix) [24, 2]. In reality, species interactions often emerge in higher-order combinations, where the relationship between two species is influenced by one or more additional species [5]. Although the importance of high-order interactions has been recognized, their impact on the stability of ecological systems has not been fully understood. Existing findings are preliminary, most of which focus on numerical simulations [20, 28].

We are particularly interested in the higher-order generalized Lotka–Volterra (HOGLV) model, designed to represent higher-order interactions within intricate ecological networks. Computing the Jacobian matrices for such high-order systems is challenging due to the curse of dimensionality, i.e., the size of variables increases exponentially with the dimensionality. To address this challenge, we propose a framework, leveraging from various tensor decompositions, including higher-order singular value decomposition (HOSVD), Canonical Polyadic decomposition (CPD), and tensor train decomposition (TTD), to improve both memory and computational efficiency in the computation of Jacobian matrix for the HOGLV model.

The article is organized as follows. In Section II, we begin with preliminaries on basic tensors. In Section III, we introduce the HOGLV model and represent it in the HOSVD, CPD, and TTD forms with memory complexity analyses. In Section IV, we derive the Jacobian matrix of the HOGLV model and offer a computational complexity analysis for each tensor decomposition-based representation. Finally, we conclude our work with discussions on future directions in Section VI.

II Tensor Preliminaries

The order of a tensor is defined as the number of its dimensions, with each dimension referred to as a mode. A $p$th-order tensor can be denoted by $\mathcal{T}\in\mathbb{R}^{n_{1}\times n_{2}\times\dots\times n_{p}}$. The tensor vector multiplication $\mathcal{T}\times_{q}\textbf{v}$ along mode $q$ for a vector $\textbf{v}\in\mathbb{R}^{n_{q}}$ is defined as $[\mathcal{T}\times_{q}\textbf{v}]_{i_{1}i_{2}\dots i_{q-1}i_{q+1}\dots i_{p}}:=\sum_{i_{q}=1}^{n_{q}}\mathcal{T}_{i_{1}i_{2}\dots i_{q}\dots i_{p}}\textbf{v}_{i_{q}}.$ The tensor matrix multiplication $\mathcal{T}\times_{q}\textbf{M}$ along mode $q$ for a matrix $\textbf{M}\in\mathbb{R}^{m\times n_{q}}$ is defined as $[\mathcal{T}\times_{q}\textbf{M}]_{i_{1}i_{2}\dots i_{q-1}ji_{q+1}\dots i_{p}}:=\sum_{i_{q}=1}^{n_{q}}\mathcal{T}_{i_{1}i_{2}\dots i_{q}\dots i_{p}}\textbf{M}_{ji_{q}}.$

In the following, we introduce three key tensor decompositions utilized in this article:

• The higher-order singular value decomposition (HOSVD) [15] of $\mathcal{T}\in\mathbb{R}^{n_{1}\times n_{2}\times\dots\times n_{p}}$ is defined as

  $$\mathcal{T}\approx\mathcal{S}\times_{1}\textbf{U}_{1}\times_{2}\textbf{U}_{2}\times_{3}\cdots\times_{p}\textbf{U}_{p},$$ (2)

  where $\mathcal{S}\in\mathbb{R}^{r_{1}\times r_{2}\times\dots\times r_{p}}$ represents the core tensor, and $\textbf{U}_{q}\in\mathbb{R}^{n_{q}\times r_{p}}$ are the factor matrices containing orthonormal columns for $q=1,2,\dots,p$. For simplicity, we may rewrite (2) in a more compact form, i.e., $\mathcal{T}\approx\mathcal{S}\times_{12\cdots p}\{\textbf{U}_{1},\textbf{U}_{2},\dots,\textbf{U}_{p}\}.$ The computation of HOSVD is numerically stable and can offer a quasi-optimal approximation.

• The Canonical Polyadic decomposition (CPD) [22] of $\mathcal{T}\in\mathbb{R}^{n_{1}\times n_{2}\times\dots\times n_{p}}$ is defined as

  $$\mathcal{T}\approx\sum_{i=1}^{r}[\textbf{T}_{1}]_{:i}\circ[\textbf{T}_{2}]_{:i}\circ\cdots\circ[\textbf{T}_{p}]_{:i},$$ (3)

  where $\textbf{T}_{q}\in\mathbb{R}^{n_{q}\times r}$ are the factor matrices for $q=1,2,\dots,p$, and $r$ is called the CP rank of $\mathcal{T}$ if it is the minimum integer achieving the decomposition. We use $\circ$ to denote the outer product. This is because CPD is not numerically stable. Moreover, the low-rank approximation of CPD is ill-posed.

• The tensor train decomposition (TTD) [26] of $\mathcal{T}\in\mathbb{R}^{n_{1}\times n_{2}\times\dots\times n_{p}}$ is defined as

  $$\mathcal{T}\approx\sum_{i_{0}=1}^{r_{0}}\cdots\sum_{i_{p}=1}^{r_{p}}[\mathcal{T}_{1}]_{i_{0}:i_{1}}\circ[\mathcal{T}_{2}]_{i_{1}:i_{2}}\circ\cdots\circ[\mathcal{T}_{p}]_{i_{p-1}:i_{p}},$$ (4)

  where $\mathcal{T}_{q}\in\mathbb{R}^{r_{q-1}\times n_{q}\times r_{q}}$ are the third-order core tensors for $q=1,2,\dots,p$, and $\{r_{0},r_{1},\dots,r_{p}\}$ is the set of TT-ranks with $r_{0}=r_{p}=1$. Similar to HOSVD, TTD is numerically stable and its low-rank approximation is quasi-optimal. Importantly, TTD is more efficient than HOSVD in both computation and memory.

III The HOGLV Model

The dynamics of complex ecological networks with higher-order interactions are often characterized by the higher-order generalized Lotka–Volterra (HOGLV) model [5, 1, 28], i.e.,

$$\begin{split}\dot{x}_{i}&=x_{i}\Big{[}r_{i}+\sum_{j=1}^{n}\textbf{A}_{ij}x_{j}+\sum_{j=1}^{n}\sum_{k=1}^{n}\mathcal{B}_{ijk}x_{j}x_{k}\\ &+\sum_{j=1}^{n}\sum_{k=1}^{n}\sum_{l=1}^{n}\mathcal{C}_{ijkl}x_{j}x_{k}x_{l}+\cdots\Big{]}\end{split}$$ (5)

for $i=1,2,\dots,n$, where $r_{i}$ is the intrinsic growth rate for species $i$, $\textbf{A}\in\mathbb{R}^{n\times n}$ is the interaction matrix whose off-diagonal elements $\textbf{A}_{ij}$ represent the effect of species $j$ upon species $i$, $\mathcal{B}\in\mathbb{R}^{n\times n\times n}$ is the third-order interaction tensor whose off-diagonal elements represent the effect that species $j$ and $k$ has upon species $i$, and $\mathcal{C}\in\mathbb{R}^{n\times n\times n\times n}$ is the fourth-order interaction tensor whose off-diagonal elements represent the effect that species $j$, $k$, $l$ has upon species $i$. In fact, the HOGLV model belongs to the family of polynomial dynamical systems.

Figure 1: High-order network representation of an ecological system with seven species, where interactions can occur between pairs or groups of nodes, and can be positive or negative.

Alternatively, we may rewrite the system of the differential equations (5) into the vector form using tensor products, which is $\dot{\textbf{x}}=\textbf{x}*[\textbf{r}+\textbf{A}\textbf{x}+\mathcal{B}\times_{23}\{\textbf{x},\textbf{x}\}+\mathcal{C}\times_{234}\{\textbf{x},\textbf{x},\textbf{x}\}+\cdots].$ The total number of the model parameters in HOGLV model, whose maximum interaction order is $M$, can be estimated as $\text{MC}_{\text{full}}\sim\mathcal{O}(n^{M})$. It can be directly seen that the $M$th-order interaction tensor has at most $n^{M}$ entries. Obviously, the memory complexity scales exponentially with the maximum order of interactions $M$. For large ecological networks with higher-order interactions, it will be challenging to handle the full representation due to severe memory limitations. In the following, we represent the HOGLV model in various tensor decomposition forms.

III-A HOSVD-based Representation

We first consider applying HOSVD to the interaction tensors to derive a HOSVD-based representation of the HOGLV model.

Proposition 1: Suppose that HOSVDs of interaction tensors $\mathcal{B}\in\mathbb{R}^{n\times n\times n}$ and $\mathcal{C}\in\mathbb{R}^{n\times n\times n\times n}$ are provided as $\mathcal{B}\approx\mathcal{B}_{0}\times_{123}\{\textbf{B}_{1},\textbf{B}_{2},\textbf{B}_{3}\}\text{ and }\mathcal{C}\approx\mathcal{C}_{0}\times_{1234}\{\textbf{C}_{1},\textbf{C}_{2},\textbf{C}_{3},\textbf{C}_{4}\},$ where $\mathcal{B}_{0}\in\mathbb{R}^{b_{1}\times b_{2}\times b_{3}}$ and $\mathcal{C}_{0}\in\mathbb{R}^{c_{1}\times c_{2}\times c_{3}\times c_{4}}$ are the core tensors with associated factor matrices $\textbf{B}_{i}\in\mathbb{R}^{n\times b_{i}}$ and $\textbf{C}_{i}\in\mathbb{R}^{n\times c_{i}}$, respectively. The HOSVD-based representation of the HOGLV model can be computed as

$$\begin{split}\dot{\textbf{x}}&=\textbf{x}*\Big{[}\textbf{r}+\textbf{A}\textbf{x}+\mathcal{B}_{0}\times_{123}\{\textbf{B}_{1},\textbf{B}_{2}^{\top}\textbf{x},\textbf{B}_{3}^{\top}\textbf{x}\}\\ &+\mathcal{C}_{0}\times_{1234}\{\textbf{C}_{1},\textbf{C}_{2}^{\top}\textbf{x},\textbf{C}_{3}^{\top}\textbf{x},\textbf{C}_{4}^{\top}\textbf{x}\}+\cdots\Big{]}.\end{split}$$ (6)

Proof:

According to the properties of tensor vector/matrix products and HOSVD, it can be shown that

$$\begin{split}&\Big{(}\mathcal{B}_{0}\times_{123}\{\textbf{B}_{1},\textbf{B}_{2},\textbf{B}_{3}\}\Big{)}\times_{23}\{\textbf{x},\textbf{x}\}\\ &=\mathcal{B}_{0}\times_{123}\{\textbf{B}_{1},\textbf{B}_{2}^{\top}\textbf{x},\textbf{B}_{3}^{\top}\textbf{x}\}.\end{split}$$

Similarly, the same principle applies to the interaction tensor $\mathcal{C}$ as well. Thus, the result follows immediately. ∎

Remark 1: Suppose that the reduced dimensions of the interaction tensors are all equal to $r$, i.e., $b_{i}=c_{i}=\cdots=r$. If the maximum order of interactions is $M$, the total number of parameters in the HOSVD-based representation of the HOGLV model can be estimated as

$$\text{MC}_{\text{hosvd}}\sim\mathcal{O}(n^{2}+M^{2}nr+r^{M}).$$

While the above memory complexity increases exponentially with the maximum order of interactions $M$, it is mitigated by the fact that the reduced dimension $r$ can be small for certain structured tensors.

III-B CPD-based Representation

We explore applying CPD to the interaction tensors to obtain a CPD-based representation of the HOGLV model.

Proposition 2: Suppose that the CPDs of the interaction tensors $\mathcal{B}\in\mathbb{R}^{n\times n\times n}$ and $\mathcal{C}\in\mathbb{R}^{n\times n\times n\times n}$ are provided as $\mathcal{B}\approx\sum_{i=1}^{b}[\textbf{B}_{1}]_{:i}\circ[\textbf{B}_{2}]_{:i}\circ[\textbf{B}_{3}]_{:i}\text{ and }\mathcal{C}\approx\sum_{i=1}^{c}[\textbf{C}_{1}]_{:i}\circ[\textbf{C}_{2}]_{:i}\circ[\textbf{C}_{3}]_{:i}\circ[\textbf{C}_{4}]_{:i},$ where $\textbf{B}_{i}\in\mathbb{R}^{n\times b}$ and $\textbf{C}_{i}\in\mathbb{R}^{n\times c}$ are the factor matrices. The CPD-based representation of the HOGLV model can be computed as

$$\begin{split}\dot{\textbf{x}}&=\textbf{x}*\Big{[}\textbf{r}+\textbf{A}\textbf{x}+\sum_{i=1}^{b}\big{(}[\textbf{B}_{2}]_{:i}^{\top}\textbf{x}\big{)}\big{(}[\textbf{B}_{3}]_{:i}^{\top}\textbf{x}\big{)}[\textbf{B}_{1}]_{:i}\\ +&\sum_{i=1}^{c}\big{(}[\textbf{C}_{2}]_{:i}^{\top}\textbf{x}\big{)}\big{(}[\textbf{C}_{3}]_{:i}^{\top}\textbf{x}\big{)}\big{(}[\textbf{C}_{4}]_{:i}^{\top}\textbf{x}\big{)}[\textbf{C}_{1}]_{:i}+\cdots\Big{]}.\end{split}$$ (7)

Proof:

From the properties of tensor vector/matrix products and CPD, it can be shown that

$$\begin{split}&\Big{(}\sum_{i=1}^{b}[\textbf{B}_{1}]_{:i}\circ[\textbf{B}_{2}]_{:i}\circ[\textbf{B}_{3}]_{:i}\Big{)}\times_{23}\{\textbf{x},\textbf{x}\}\\ &=\sum_{i=1}^{b}\big{(}[\textbf{B}_{2}]_{:i}^{\top}\textbf{x}\big{)}\big{(}[\textbf{B}_{3}]_{:i}^{\top}\textbf{x}\big{)}[\textbf{B}_{1}]_{:i}.\end{split}$$

Likewise, the same principle applies to the interaction tensor $\mathcal{C}$ as well. Therefore, the result follows immediately. ∎

Remark 2: Suppose that the CP ranks of the interaction tensors are all equal to $r$, i.e., $b=c=\dots=r$. If the maximum order of interactions is $M$, the total number of parameters in the CPD-based representation of the HOGLV model can be estimated as

$$\text{MC}_{\text{cpd}}\sim\mathcal{O}(n^{2}+M^{2}nr).$$

Remarkably, the above memory complexity increases quadratically with both the maximum order of interactions $M$ and the number of species $n$. Furthermore, when the CP rank $r$ is small, the CPD-based representation of the HOGLV model can significantly enhance memory efficiency. However, it is important to note that CPD may not be numerically stable, which may limit the applicability of this representation.

III-C TTD-based Representation

We investigate the use of TTD on the interaction tensors to derive a TTD-based representation of the HOGLV model.

Proposition 3: Suppose that the TTDs of the interaction tensors $\mathcal{B}\in\mathbb{R}^{n\times n\times n}$ and $\mathcal{C}\in\mathbb{R}^{n\times n\times n\times n}$ are provided as $\mathcal{B}\approx\sum_{i_{0}=1}^{b_{0}}\cdots\sum_{i_{3}=1}^{b_{3}}[\mathcal{B}_{1}]_{i_{0}:i_{1}}\circ[\mathcal{B}_{2}]_{i_{1}:i_{2}}\circ[\mathcal{B}_{3}]_{i_{2}:i_{3}}\text{ and }\mathcal{C}\approx\sum_{i_{0}=1}^{c_{0}}\cdots\sum_{i_{4}=1}^{c_{4}}[\mathcal{C}_{1}]_{i_{0}:i_{1}}\circ[\mathcal{C}_{2}]_{i_{1}:i_{2}}\circ[\mathcal{C}_{3}]_{i_{2}:i_{3}}\circ[\mathcal{C}_{4}]_{i_{3}:i_{4}},$ where $\mathcal{B}_{i}\in\mathbb{R}^{b_{i-1}\times n\times b_{i}}$ and $\mathcal{C}_{i}\in\mathbb{R}^{c_{i-1}\times n\times c_{i}}$ are the third-order core tensors. The TTD-based representation of the HOGLV model can be computed as

$$\small\begin{split}\dot{\textbf{x}}&=\textbf{x}*\Big{[}\textbf{r}+\textbf{A}\textbf{x}+\sum_{i_{0}=1}^{b_{0}}\cdots\sum_{i_{3}=1}^{b_{3}}([\mathcal{B}_{2}]_{i_{1}:i_{2}}^{\top}\textbf{x})([\mathcal{B}_{3}]_{i_{2}:i_{3}}^{\top}\textbf{x})[\mathcal{B}_{1}]_{i_{0}:i_{1}}\\ +&\sum_{i_{0}=1}^{c_{0}}\cdots\sum_{i_{4}=1}^{c_{4}}([\mathcal{C}_{2}]_{i_{1}:i_{2}}^{\top}\textbf{x})([\mathcal{C}_{3}]_{i_{2}:i_{3}}^{\top}\textbf{x})([\mathcal{C}_{4}]_{i_{3}:i_{4}}^{\top}\textbf{x})[\mathcal{C}_{1}]_{i_{0}:i_{1}}+\cdots\Big{]}.\end{split}$$ (8)

Proof:

The proof is similar to Proposition 2. ∎

Remark 3: Suppose that the TT-ranks of the interaction tensors are all equal to $r$, i.e., $b_{i}=c_{i}=\cdots=r$. Note that the first and last TT-ranks of each tensor are always equal to 1. If the maximum order of interactions is $M$, the total number of parameters in the TTD-based representation of the HOGLV model can be estimated as

$$\text{MC}_{\text{ttd}}\sim\mathcal{O}(n^{2}+M^{2}nr^{2}).$$

While the above memory complexity is slightly higher than that of the CPD-based representation, it remains significantly more efficient than the full and HOSVD-based representations. Importantly, TTD can offer the crucial advantages of numerical stability and quasi-optimal approximation.

IV Complexity Analyses of Linear Stability

Linear stability analysis has proven to be a successful tool in the study of complex ecological networks, enhancing our understanding of the intricate relationships between stability and biodiversity [3, 4, 24]. A plethora of research has focused on the linear stability of the GLV model [19, 6]. The essence of linear stability lies in the computation of the Jacobian matrix. Here, we derive the Jacobian matrix for the HOGLV model and perform a computational complexity analysis of computing the Jacobian matrix using the full, HOSVD-based, CPD-based, and TTD-based representations.

Proposition 4: Suppose that the equilibrium point of the HOGLV model is $\textbf{x}^{*}$. The Jacobian matrix of the HOGLV model evaluated at $\textbf{x}^{*}$ can be computed as

$$\begin{split}\textbf{M}&=\textbf{X}^{*}\Big{[}\textbf{A}+\mathcal{B}\times_{2}\textbf{x}^{*}+\mathcal{B}\times_{3}\textbf{x}^{*}+\mathcal{C}\times_{23}\{\textbf{x}^{*},\textbf{x}^{*}\}\\ &+\mathcal{C}\times_{24}\{\textbf{x}^{*},\textbf{x}^{*}\}+\mathcal{C}\times_{34}\{\textbf{x}^{*},\textbf{x}^{*}\}+\cdots\Big{]},\end{split}$$ (9)

where $\textbf{X}^{*}\in\mathbb{R}^{n\times n}$ is a diagonal matrix that contains $\textbf{x}^{*}$ along its diagonal.

Proof:

Based on the properties of tensor vector products, it can be demonstrated that

$$\frac{d}{d\textbf{x}}\mathcal{B}\times_{23}\{\textbf{x},\textbf{x}\}=\mathcal{B}\times_{2}\textbf{x}+\mathcal{B}\times_{3}\textbf{x}.$$

Likewise, the same derivative principle is applicable to the interaction tensor $\mathcal{C}$. Consequently, the result can be derived in a manner similar to the computation of the Jacobian matrix for the GLV model [8]. ∎

Remark 4: If the maximum order of interactions is $M$, the computational complexity of computing the Jacobian matrix of the HOGLV model can be estimated as

$$\text{CC}_{\text{full}}\sim\mathcal{O}(M^{2}n^{M}).$$

Similar to its memory complexity, the computational complexity of the full representation scales exponentially with the maximum order of interactions $M$. Therefore, computing the Jacobian matrix for large ecological networks with higher-order interactions is an exceptionally challenging task. To address this challenge, we can utilize tensor decomposition-based representations of the HOGLV model, which significantly reduces the computational complexity. In the following, we estimate the computational complexity of computing the Jacobian matrix using the HOSVD-based, CPD-based, and TTD-based representations.

Remark 5: Given the HOSVD-based representation of the HOGLV model (6) with reduced dimension $r$ for all interaction tensors, if the maximum order of interactions is $M$, the computational complexity of computing the Jacobian matrix can be estimated as

$$\text{CC}_{\text{hosvd}}\sim\mathcal{O}(M^{2}n^{2}r+M^{3}nr+M^{2}nr^{M}).$$

Remark 6: Given the CPD-based representation of the HOGLV model (7) with CP rank $r$ for all interaction tensors, if the maximum order of interactions is $M$, the computational complexity of computing the Jacobian matrix can be estimated as

$$\text{CC}_{\text{cpd}}\sim\mathcal{O}(M^{2}n^{2}r+M^{3}nr).$$

Remark 7: Given the TTD-based representation of the HOGLV model (8) with TT-rank $r$ for all interaction tensors (except for the first and last TT-ranks), if the maximum order of interactions is $M$, the computational complexity of computing the Jacobian matrix can be estimated as

$$\text{CC}_{\text{ttd}}\sim\mathcal{O}(Mnr^{M}+n^{2}r^{M}).$$

The proofs of the above remarks can be found in the appendix. The computational complexity of computing the Jacobian matrix using the HOSVD-based, CPD-based, and TTD-based representations exhibits a similar behavior to their memory complexity. Particularly, the CPD-based representation stands out with the lowest computational complexity, which does not grow exponentially with the maximum order of interactions $M$. Nonetheless, accurately computing CPD becomes challenging for tensors with larger dimensions. Conversely, the TTD-based representation, although having a higher computational complexity compared to the CPD-based representation, outperforms the full and HOSVD-based representations with guaranteed numerical stability.

After obtaining the Jacobian matrix, we can compute its eigenvalues to determine the linear stability of the HOGLV model at an equilibrium point. It is important to acknowledge that computing the eigenvalues of the Jacobian matrix can be computationally demanding for large ecological networks, but this specific aspect is not the primary focus of this article. Nevertheless, recent research [23, 16] has explored the use of TTD to expedite eigenvalue computations, which can be useful here for future investigations.

V Numerical Examples

We present three case studies to illustrate our framework. All three examples were conducted on a Macintosh machine equipped with 32 GB RAM and an Apple M1 Pro chip (10-core CPU at 3.2 GHz, 16-core GPU, and 16-core Neural Engine) using MATLAB R2022a with the Tensor Toolbox 3.0 [21] and TT Toolbox [25].

Figure 2: Two ecological networks, fully characterized by their pairwise and third-order interactions, with positive interactions colored in blue and negative interactions in red.

V-A Computing Jacobian Matrices

We computed the Jacobian matrices of two ecological networks captured by the HOGLV model using (9). For simplicity, we set the total number of species and the maximum order of interactions to 4 and 3, respectively. The structures of the two ecological networks are presented in Fig. 2, which share pairwise interactions. Suppose that the equilibrium point is $\textbf{x}^{*}=\begin{bmatrix}1&1&1&1\end{bmatrix}^{\top}$ for both networks. The associated Jacobian matrices can be computed as

$$\textbf{M}_{1} =\begin{bmatrix}-0.1000&-0.2147&\phantom{-}1.0826&\phantom{-}0.8781\\ \phantom{-}0.1468&-0.1000&-0.6852&\phantom{-}0\\ \phantom{-}0&\phantom{-}0.0923&-0.1000&\phantom{-}0\\ \phantom{-}0&\phantom{-}0&\phantom{-}0&-0.1000\end{bmatrix},$$

$$\textbf{M}_{2} =\begin{bmatrix}-0.1000&-0.6236&-1.0826&-0.8781\\ \phantom{-}0.1468&-0.1000&-0.6852&\phantom{-}0\\ \phantom{-}0&\phantom{-}0.0923&-0.1000&\phantom{-}0\\ \phantom{-}0&\phantom{-}0&\phantom{-}0&-0.1000\end{bmatrix}.$$

The first Jacobian matrix $\textbf{M}_{1}$ has one positive eigenvalue, indicating that the first ecological network is locally unstable around the equilibrium point. The second Jacobian matrix $\textbf{M}_{2}$ has all eigenvalues negative, suggesting that the second ecological network is locally stable near the equilibrium point. This example demonstrates that higher-order interactions can significantly influence the stability of complex ecological networks.

V-B Memory Comparison

We conducted a comparison of the total number of model parameters for the full, HOSVD-based, CPD-based, and TTD-based representations of the HOGLV model using random sparse tensors with varying dimensions. We set the maximum order of interactions to 4 and randomly generated sparse interaction tensors $\mathcal{B}\in\mathbb{R}^{n\times n\times n}$ and $\mathcal{C}\in\mathbb{R}^{n\times n\times n\times n}$ for $n=10,20,\dots,50$. By employing the HOSVD, CPD, and TTD on the interaction tensors, we computed their corresponding representations of the HOGLV model. The results are presented in Table I. The HOSVD-based representation is effective in reducing the number of model parameters for relatively small dimensions ($n=10,20$). However, as the dimension increases, its number of model parameters surpasses that of the full representation. Moreover, the CPD-based representation demonstrates the lowest memory usage. However, we failed to accurately compute the CPDs of the interaction tensors for larger dimensions ($n=30,40,50$) due to the issue of numerical instability. Thus, the TTD-based representation can be considered as a robust and viable choice, which maintains a reasonably low number of model parameters compared to the full representation.

TABLE I: Comparison of the total number of model parameters for the full, HOSVD-based, CPD-based, and TTD-based representations of the HOGLV model across varying system dimensions.

Dimension	$10$	$20$	$30$	$40$	$50$
Full	11,110	168,420	837,930	2,625,640	6,377,550
HOSVD	810	131,608	844,230	2,636,840	6,395,050
CPD	460	6,860	-	-	-
TTD	890	36,880	238,530	920,040	2,687,550

V-C Computation Comparison

We compared the computational time required for computing the Jacobian matrix of the HOGLV model using the full representation and the TTD-based representation. We exclusively considered the TTD-based representation here due to its numerical stability and low computational complexity. More importantly, it benefits from a well-established TT-algebra, making it a suitable choice for this computation. We set the maximum order of interactions to 4 and randomly generated interaction tensors $\mathcal{B}\in\mathbb{R}^{n\times n\times n}$ and $\mathcal{C}\in\mathbb{R}^{n\times n\times n\times n}$ in the TTD format with low TT-ranks for varying dimensions. In other words, we assumed that the HOGVL model is initially provided in the TTD-based representation. The findings are displayed in Fig. 3. Evidently, the TTD-based representation offers a substantial time advantage in computing the Jacobian matrix compared to the full representation. It is noteworthy that the full representation encounters memory limitations for $n\geq 275$, which fails to complete the computation. On the other hand, the TTD-based representation maintains computational efficiency even for larger dimensions, as demonstrated in Table II.

TABLE II: Computation time for calculating the Jacobian matrix of the TTD-based representation of the HOGLV model for large system dimensions, where the full representation fails due to memory constraints.

Dimension	$1000$	$2000$	$3000$	$4000$	$5000$
Time (s)	0.0259	0.1675	0.5304	1.1228	1.9990

Figure 3: Comparison of computational time for computing the Jacobian matrix between the full and TTD-based representations.

VI Conclusion

In this article, we proposed a framework to tackle the problem of complex ecological systems with high-order interactions. In particular, we were interested in determining the linear stability of the HOGLV model by computing the Jacobian matrix and its eigenvalues. A crucial challenge in the analysis of complex systems and coupled representations is their exponentially increasing size. In our approach, we explored various tensor decomposition techniques for good scalability with detailed memory and computational complexity investigations.

The stability analysis of dynamical systems with high-order interactions is not limited to ecological systems, but also appears in various fields, such as opinion dynamics and social networks, where the aggregated impact of high-order interactions plays a significant role. Another area of interest is the estimation of interaction matrices and tensors of different orders [17]. Additionally, the study of high-order time-varying systems, where the interaction tenors are time-varying, is a promising direction for future research.

We provide detailed proofs for Remarks 4 and 5 here. Remarks 6 and 7 can be proved in a similar manner.

Proof of Remark 4

We first consider the computational complexity of tensor vector multiplications, which can be estimated as

$$\mathcal{O}(2n^{3}+3n^{4}+\cdots\binom{M}{2}n^{M})\sim\mathcal{O}(M^{2}n^{M}).$$

Additionally, the computational complexity of matrix additions and multiplications is much less than $\mathcal{O}(M^{2}n^{M})$. Therefore, the overall complexity of computing the Jacobian matrix of the HOGLV model using the full representation is about $\mathcal{O}(M^{2}n^{M})$.

Proof of Remark 5

Suppose that the reduced dimensions of the interaction tensors are all equal to $r$, i.e., $b_{i}=c_{i}=\cdots=r$. We first estimate the computational complexity of computing the term $\mathcal{B}_{0}\times_{123}\{\textbf{B}_{1},\textbf{B}_{2}^{\top}\textbf{x},\textbf{B}_{3}\}$, which is about

$$\mathcal{O}(nr+nr^{3}+nr^{2}+n^{2}r).$$

Therefore, the computational complexity of tensor vector multiplications can be approximated by

$$2\mathcal{O}(nr+nr^{3}+nr^{2}+n^{2}r)$$

$$+ 3\mathcal{O}(2nr+nr^{2}+nr^{3}+nr^{4}+n^{2}r)$$

$$+ \cdots$$

$$+ \binom{M}{2}\mathcal{O}((M-1)nr+nr^{2}+\cdots+nr^{M}+n^{2}r)$$

$$\sim\mathcal{O}(M^{3}nr+M^{2}nr^{M}+M^{2}n^{2}r).$$

Again, the computational complexity of matrix additions and multiplications can be negligible. Thus, the overall complexity of computing the Jacobian matrix of the HOGLV model using the HOSVD representation is about $\mathcal{O}(M^{3}nr+M^{2}nr^{M}+M^{2}n^{2}r)$.