WGMwgm QEqe EPep PMSpms BECbec DEde

[1]

[1]This document is the results of the research project funded by Non-
commercial Foundation for the Advancement of Science and Education INTELLECT.

[type=editor, auid=000,bioid=1, orcid=0000-0002-0513-4425]

[1]

[1]

Conceptualization of this study, Methodology, Software, Data curation, Writing
- Original draft preparation

1]organization=HSE University, addressline=20 Myasnitskaya St., city=Moscow,
postcode=101000, country=Russia

2]organization=Institute of Astronomy RAS, addressline=48 Pyatnitskaya St.,
city=Moscow, postcode=119017, country=Russia

[cor1]Corresponding author

# Machine learning methods for the search for L&T brown dwarfs in the data of
modern sky surveys.

Aleksandra Avdeeva<EMAIL_ADDRESS>[  [

###### Abstract

According to various estimates, brown dwarfs (BD) should account for up to 25
percent of all objects in the Galaxy. However, few of them are discovered and
well-studied, both individually and as a population. Homogeneous and complete
samples of brown dwarfs are needed for these kinds of studies. Due to their
weakness, spectral studies of brown dwarfs are rather laborious. For this
reason, creating a significant reliable sample of brown dwarfs, confirmed by
spectroscopic observations, seems unattainable at the moment. Numerous
attempts have been made to search for and create a set of brown dwarfs using
their colours as a decision rule applied to a vast amount of survey data. In
this work, we use machine learning methods such as Random Forest Classifier,
XGBoost, SVM Classifier and TabNet on PanStarrs DR1, 2MASS and WISE data to
distinguish L and T brown dwarfs from objects of other spectral and luminosity
classes. The explanation of the models is discussed. We also compare our
models with classical decision rules, proving their efficiency and relevance.

###### keywords:

brown dwarfs machine learning surveys

##  1 Introduction

Brown dwarfs are substellar objects that were theoretically predicted (Kumar,
1963; Hayashi and Nakano, 1963) and then discovered 30 years later (Rebolo et
al., 1995; Nakajima et al., 1995). Since then, the search (Luhman, 2013;
Burningham et al., 2013; Carnero Rosell et al., 2019) and systematic study of
known brown dwarfs (Kirkpatrick et al., 1999; Skrzypek et al., 2016;
Kirkpatrick et al., 2021) has not stopped. Their mass is insufficient to start
and maintain stable hydrogen fusion, which causes them to cool over time. The
peak of the radiation intensity falls into the infrared range, so objects are
rather faint in the visible spectrum. In the spectral classification, brown
dwarfs occupy spectral types L, T and Y.

According to studies (Mužić et al., 2017), the number of brown dwarfs in the
Galaxy ranges from 25 to 100 billion objects (with the total number of objects
ranging from 100 to 500 billion). Homogeneous and complete samples of brown
dwarfs are needed for various kinds of studies: kinematic studies (Smith et
al., 2014), studies of binary stars with brown dwarfs (Lodieu et al., 2014),
and studies of the parameters of the Galaxy. Brown dwarfs occupy the boundary
between stars and planets, and studying their properties helps to refine our
understanding of this boundary. Complete and uniform catalogues enable to
identify and characterize brown dwarfs with greater accuracy, allowing for a
better determination of the lower mass limit for stellar formation and the
upper mass limit for planet formation. Moreover, brown dwarfs share
similarities with giant exoplanets, making them valuable analogs for studying
exoplanetary atmospheres. By studying the atmospheres of brown dwarfs, similar
to exoplanets, we can gain insights into the processes and conditions that
govern exoplanet atmospheres, including the presence of clouds, atmospheric
composition, and thermal profiles.

Probably the most topical issue regarding brown dwarfs is the L/T transition
(Artigau et al., 2009; Gillon et al., 2013; Khandrika et al., 2013). The L/T
transition in brown dwarfs is a fascinating phenomenon that is characterized
by a sharp change in the near-infrared colours and brightness of brown dwarfs.
It is believed to be driven by several possible mechanisms. Cloud models link
the sharp transition to the sinking of dust clouds below the photosphere
(Marley and Robinson, 2015; Charnay et al., 2018). Instability in the carbon
chemistry of brown dwarf atmospheres has been proposed as another mechanism
contributing to the L/T transition (Tremblin et al., 2019). Adiabatic
convection triggered by this instability can lead to variability across the
L/T spectral sequence. Cloud dispersal has emerged as a potential mechanism
for the L/T transition (Tan and Showman, 2019). It has been suggested that
clouds with larger particle sizes dissipate more easily than clouds with
smaller particle sizes. The shift from L to T spectral types may be
accompanied by a change from small to large particles, leading to the
fragmentation of clouds and a transition to cloud-free atmospheres (Burningham
et al., 2017; Saumon and Marley, 2008). A detailed overview of the problem can
be found in Vos et al. (2019).

To uncover the nature of the phenomenon mentioned above, as well as to refine
the statistical characteristics of brown dwarfs, complete and uniform
catalogues of brown dwarfs are required. While spectroscopy is essential for
confirming the nature of a brown dwarf and studying its detailed properties,
conducting spectroscopic observations for a large number of objects across the
entire sky is time-consuming and resource-intensive. Photometric surveys, on
the other hand, can cover a much larger area of the sky and efficiently
capture data on numerous celestial objects simultaneously.

By employing colour selection techniques in photometric surveys, one can
identify objects that exhibit colours indicative of brown dwarf
characteristics. The advantage of using photometric surveys is that they allow
for systematic and wide-scale screening of potential brown dwarf candidates,
helping to identify promising targets for subsequent spectroscopic
observations.

As an illustration, Skrzypek et al. (2016) accomplished this feat by
effectively employing data from three surveys: SDSS, UKIDSS, and WISE. They
employed a specific colour selection criterion, namely $(Y-J)_{Vega}>0.8$ and
$J<17.5$, as a decision rule. Through this approach, they were able to discern
approximately 1300 brown dwarfs within an area of 3000 square degrees, which
corresponds to approximately 7.5% of the celestial sphere. An additional
noteworthy implementation in the quest for brown dwarfs is exemplified by
Carnero Rosell et al. (2019). Their study also incorporated a decision rule
utilizing data from the DES, VHS, and WISE surveys. The following criteria
were applied: $(i-z)>1.2$, $(z-Y)>0.15$, $(Y_{AB}$-$J_{Vega})>1.6$, and
$z<22$. The imposition of a magnitude limit on the $z$ band was necessary to
ensure the completeness of the dataset, thereby excluding any missing values.
Within an area spanning 2400 square degrees about 5.8% of the celestial
sphere, approximately 12 thousand brown dwarfs were successfully identified
through their approach. Carnero Rosell et al. (2019) also presents a
comprehensive review of other colour selection works.

When it comes to using colour selection to search for brown dwarfs,
incorporating machine learning methods can provide significant advantages.
Machine learning techniques can enhance the effectiveness and efficiency of
the colour selection process by leveraging large datasets and complex
algorithms to identify patterns and make more accurate predictions.

Machine learning methods can help uncover subtle relationships and
correlations in multi-dimensional colour space, allowing for the
identification of distinct colour signatures associated with brown dwarfs.
This can be particularly valuable when dealing with complex and overlapping
colour distributions between different objects.

ML methods are increasingly being used for classifying astronomical objects
due to the vast amount of data collected in the past decades. For example,
Maravelias et al. (2022) combined Support Vector Machine (SVM), Random Forest
(RF) and Multilayer Perceptron to classify massive stars in nearby galaxies.
The accuracy of the test dataset was 83%. Applying it to other galaxies (not
included in the dataset), namely IC 1613, WLM and Sextans A, achieved an
accuracy result at the level of 70%, which the authors attribute to different
metallicity and extinction effects. The missing data were filled in with
simple averages and the Iterative Imputer method of the Scikit-learn library
(Pedregosa et al., 2011). The Iterative Imputer method calculates the missing
values based on the present values for the features in the same manner as
regression models do. This method, being at the same time more robust, showed
a better performance in the work.

The interpretable machine learning techniques (Localized General Matrix LVQ
and RF) were used in Mohammadi et al. (2022) to detect extragalactic Ultra-
compact dwarfs and Globular Clusters. Authors analysed the importance of
features and compared them with features that carry physical information of
the objects.

This work aims to develop an additional tool for the search for brown dwarfs
in large photometric surveys with machine learning methods. That is, based on
the set of magnitudes and colours of the object, the model must determine
whether the given object is a brown dwarf or not. We also compare our results
with some classical decision rules: Burningham et al. (2013) and Carnero
Rosell et al. (2019). The summary of these rules is shown in Tab. 2. The tool
will be used in future work for the search for previously undiscovered brown
dwarfs.

This paper is organised as follows. Sec. 2 describes the dataset and
preprocessing of the data, including feature engineering, augmentation and the
approach to handling the missing values. In Sec. 3 we apply the machine
learning methods to the dataset and explore the importance of features.
Finally, in Sec. 4 we compare the performance of machine learning models and
the classical decision rules and discuss the robustness of the models.

##  2 Building the dataset and preprocessing

The dataset is based on L and T brown dwarfs from the Best et al. (2018)
catalogue. The catalogue contains information on 1601 L and T type brown
dwarfs and 8287 M type red dwarfs, the spectral class closest in physical
characteristics to brown dwarfs. Magnitudes in 12 photometric bands and their
errors are provided: $g,r,i,z,y$ of Pan-STARRS 1 (Chambers et al., 2016),
$J,H,Ks$ of 2MASS space mission (Cutri et al., 2003) and $W1,W2,W3,W4$ of WISE
space mission (Cutri et al., 2021).

Figure 1: Response curves of photometric systems.

The catalogue also contains astrometrical information: position, parallax and
proper motion. In addition, references to the literature are presented, from
which data on proper motion and parallax was taken.

For our Machine Learning models we consider brown dwarfs to be a positive
class. To create a representative distribution of negative class objects, we
examined the distribution of 100 thousand stars from Gaia DR3 (Collaboration
et al., 2016, 2022) by absolute magnitude $M_{G}$ (Fig. 2a). We have selected
1791 objects from A0 to K9 spectral class in the proportions observed in Fig.
2a from the Simbad111http://simbad.cds.unistra.fr/simbad/ database of
astronomical objects. Objects that are presented in Simbad are usually well-
studied and have solid spectral classifications. Gaia data seem to be short
regarding M-type dwarfs, especially, after M3, so we adopt their distribution
from Best et al. (2018). The distribution of the dataset we have obtained is
shown in Fig. 2b.

It should be noted, the obtained distribution is most likely incomplete, both
in M-dwarfs and in earlier spectral classes. The deficiency in M-dwarfs could
be due to Best et al. (2018) being not complete for M-type dwarfs outside of
10-pc. The objects of A0 to K9 might tend to be brighter than the actual
population and in fact have an underdensity in the range between $12^{m}$ and
$15^{m}$ in $i_{PS1}$ (Fig 3a), which we, however, assume to be not relevant
in case of using only colour indices as features.

Figure 2: Absolute colour distribution for the Gaia sample (a) and the dataset
used in this work (b).

The objects selected from Simbad were cross-matched with data from the Pan-
STARRS DR1, 2MASS, and ALLWISE catalogues. We chose the matching radius to be
$1^{\prime\prime}$, which is a reasonable value for most surveys for objects
with low proper motions, including those used in this paper.

The resulting dataset contains 5669 objects, 1601 of which are of the positive
target class. The dataset is available
online222https://github.com/iamaleksandra/ML-Brown-Dwarfs/. As the peak
intensity of brown dwarfs falls in the infrared part of the spectrum, their
magnitudes in the optical photometric bands (bands $g$,$r$,$i$) are most
likely almost beyond the sensitivity limit of the telescope and therefore are
missing from the data. The $g$ and $r$ bands of the Pan-STARRS data are
missing for almost all objects, so we don’t use these magnitudes as features.
The $i$ band values are missing for about one-third of the objects of the
positive class and a small number of objects of the negative class. Magnitude
values in this band are important to us, also for comparison with the
classical rules, so we keep them. We also remove $W3$ and $W4$ as they are of
poor quality, 90th percentile of the error of magnitude measurement for both
magnitudes is about half of a magnitude, which is very noisy for the
classification problem we want to solve. As a result, we have 7 magnitudes for
each object: $i_{PS1}$, $z_{PS1}$, $y_{PS1}$, $J$, $H$, $Ks$, $W1$, $W2$.

###  2.1 Train, validation and test

The data is divided into training (60%), validation (20%) and test (20%) sets
in a stratified fashion using using sci-kit learn train_test_split method. The
features are scaled using StandardScaler and model hyperparameters are
selected on the validation set using optuna (Akiba et al., 2019). The final
performance of the model is calculated on the test set.

###  2.2 Augmentation

As the number of objects of negative class is almost 2.5 times higher than the
number of positive class objects, we use augmentation with Gaussian-
distributed noise as an oversampling method to make the dataset more balanced.
The data was divided into training, validation and test subsamples prior to
augmentation in order to not let models see the augmented data in the test
sample, while being trained on the original prototype of these augmented
objects.

We augment the data of positive and negative classes separately, using all
objects of positive class and only objects in the range between $12^{m}$ and
$15^{m}$ in $i_{PS1}$ in negative class objects. For each feature error, we
calculate the mean and the standard deviation. We then generate normally
distributed noise values with the same parameters of the distribution. The
noise is added to all values of the corresponding features, which do not have
any missed values. Thus we have 8364 objects of which 4155 are positive and
4209 are negative.

###  2.3 Feature engineering

In a classification of astronomical objects, as in various types of
astronomical problems, the colours of objects are even more important than the
magnitudes. Colours are characteristic of the energy distribution in the
spectrum and are almost independent of distance. To take this into account
when classifying, we have added several features - colour indices:
$(i-z)_{PS1}$, $(i-y)_{PS1}$, $(z-y)_{PS1}$, $z_{PS1}-J$, $y_{PS1}-J$, $J-H$,
$H-Ks$, $Ks-W1$, $W1-W2$. They are also often used as the colour selection to
distinguish brown dwarfs from other objects. We only use colours of the most
spectroscopically adjacent magnitudes (see Fig. 1). The two exceptions are
$z_{PS1}-J$ since it is commonly used as a colour selection and proved to be
useful, and $(i-y)_{PS1}$ since it turned out to be extremely useful for
classification (see Sec. 4 for details).

After this procedure, the table contains 17 features for each of the 8364
objects, listed in Tab. 1. Fig. 3 shows how objects of the target class look
compared to objects of all other classes in a two-dimensional slice of feature
space. We use all of the magnitudes and colours simultaneously (even though
the latter ones are the linear combination of the first ones by design) since
we process the missing values independently, which sometimes violates these
relations (for the details see Sec. 2.4).

As one can see (Fig. 3a), the upper limit of magnitude $i_{PS1}$ differs a lot
for our positive and negative objects. This is due to the procedure of
building the dataset, i.e. Simbad has four times more objects with
$i_{SDSS}>20$ and a spectral type $\geq L0$ than objects with a spectral type
$<L0$. Catalogues, however, contain large amount of faint objects, the
majority of which are not brown dwarfs. We therefore should avoid models based
primarily on PS1 magnitudes. Although our dataset can be called balanced
regarding other magnitudes (from 2MASS and WISE), the magnitude is not only a
function of the luminosity of an object, but also a function of the distance
to it, so it is not advisable to rely on these magnitudes as well.

Thus, three cases for each approach are investigated: all magnitudes and
colours are used as features (we call this case ”all features”), no Pan-STARRS
magnitudes are used (”w/o PS magnitudes”) and no magnitudes used at all (”only
colours”).

Figure 3: Augmented data: objects of different classes on the magnitude-
magnitude (a) and colour-magnitude (b) diagrams. Brown triangles are the
objects of positive class (L&T brown dwarfs), blue squares are the objects of
negative class. Here blue points are plotted on top of the brown triangles.

###  2.4 Processing of missing values

As it was mentioned earlier, the dataset has a significant number of missing
values. At the part of the spectrum with longer wavelength, this is most
likely due to the sensitivity limit of the telescope: brown dwarfs are rather
faint objects, and their emission maximum falls into the infrared part of the
spectrum. Missing values at the shorter wavelength part of the spectrum seem
to occur due to poor-quality measurements or artefacts.

To deal with the missing values, Maravelias et al. (2022) used imputing by
means and the Iterative Imputer of the Sckit-learn library, showing that
Iterative Imputer provides results more robust and effective in terms of
classification.

We test the method by additionally throwing out the magnitude values for 5
percent of objects, which do not have missing values of the particular
feature. Then we impute these values using the Iterative Imputer and compare
the results to the original feature values for the object.

Tab.1 shows the results of the imputation with the following parameters of the
Iterative Imputer:

⬇

estimator=ExtraTreesRegressor

(n_estimators=150

max_features=14

max_depth=15

min_samples_split=12

initial_strategy=’median’

max_iter=20)

Tab. 1 contains the information about the fraction of missing values of a
particular feature in the dataset and the number of objects that were withheld
for the testing. We also compare the 90th percentile error of the feature
measurement (the error of magnitude is usually given in the catalogue and the
error of the colour index is calculated as the square root of the sum of the
squared errors of the magnitudes involved) with the 90th percentile
discrepancy in the actual value of the feature and the value predicted by the
imputer.

In most cases, the 90th percentile for the discrepancy between the imputed
values and the original ones is compatible with the 90th percentile for the
error of the feature. Even though the colour values calculated in the Sec. 2.3
are directly related to the magnitude values, it was decided to apply
Iterative Imputer to them independently. This allows one to achieve better
results and avoid large errors in the calculation of colour values, as can be
seen from Tab. 1.

In Figure 4, we can see an example of imputing. The top panel displays a
magnitude-magnitude diagram containing both the original and imputed data.
It’s worth noting that a significant portion of the missing data points are
located in the fainter range of the $i$ magnitude. The bottom panel compares
the original data to the imputed values for the same objects. Although there
are some discrepancies up to $0.5^{\textrm{m}}$ in the $y_{PS1}-J$ color, most
of the values are predicted accurately. Specifically, for 90% of the stars,
$y_{PS1}-J$ value have an error less than $0.063^{\textrm{m}}$.

Table 1: Properties of the dataset and the results of the imputing test for every feature for the train part of the dataset. The table includes the fraction of missing values in the dataset (after augmentation) and the number of objects which were withheld for the test (5% of the objects for which the value of the feature was presented). We compare the 90th percentile for the discrepancy between the imputed values and the original ones with the 90th percentile for the error in the measurement of the feature value. The error of the colour index is calculated as the square root of the sum of the squared errors of the magnitudes involved. Feature | Fraction of | Number of | 90th percentile | 90th percentile
---|---|---|---|---
| missing values | objects withheld | of the error | of discrepancy
iPS1 | 17% | 208 | 0.050 | 0.070
zPS1 | 5.5% | 237 | 0.050 | 0.091
yPS1 | 2.2% | 245 | 0.060 | 0.109
J | 8.8% | 228 | 0.100 | 0.101
H | 8.8% | 228 | 0.120 | 0.107
Ks | 8.9% | 228 | 0.110 | 0.075
W1 | 2.1% | 245 | 0.141 | 0.085
W2 | 2.0% | 242 | 0.080 | 0.096
(i-z)PS1 | 18.4% | 204 | 0.067 | 0.050
(i-y)PS1 | 18.0% | 205 | 0.072 | 0.038
(z-y)PS1 | 6.7% | 234 | 0.073 | 0.077
zPS1-J | 12.3% | 219 | 0.122 | 0.048
yPS1-J | 10.9% | 223 | 0.130 | 0.063
J-H | 8.8% | 228 | 0.164 | 0.154
H-Ks | 9.1% | 225 | 0.170 | 0.145
Ks-W1 | 9.9% | 226 | 0.183 | 0.143
W1-W2 | 2.2% | 245 | 0.168 | 0.176

Figure 4: An example of imputed values for y${}_{PS1}-$J and yPS1 features.
(Top) Original data (pink circles), initially missing data (green plus) and
manually missing values (blue crosses) on the colour-magnitude diagram.
(Bottom) Comparison of the original data (pink circles) and the values,
imputed via the Iterative Imputer method (white crosses). An index is used to
avoid contamination from the imputed values of other features.

##  3 Application of machine learning

Four approaches are tested during the work: Random Forest (RF), Support Vector
Machines (SVM), XGBoost and TabNet. As it was said in Sec. 2, we investigate
three cases for each approach: ”all features”, ”w/o PS magnitudes” and ”only
colours”. We calculate the score and feature importance for each model, using
SHAP (Lundberg and Lee, 2017) for RF, XGBoost, SVM and TabNet. Although TabNet
has built-in capabilities of calculating importance of features,based on the
attention mechanism’s dynamic selection of input features, we use SHAP as well
so as we could compare the results properly.

The SHAP method works by evaluating the model’s prediction for each instance
while permuting the values of a specific feature. This involves shuffling the
values of the chosen feature while keeping the rest of the features unchanged.
The difference between the model’s prediction with the original feature values
and the prediction with the permuted feature values is used to calculate the
Shapley value.

We compare all models to the classical decision rules Carnero Rosell et al.
(2019) and Burningham et al. (2013). The Matthews correlation coefficient
(MCC) is chosen as the primary metric since it takes into account both False
Positive and False Negative predictions. It can be calculated as follows:

$\textsc{MCC}=\frac{\textsc{TP}\times\textsc{TN}-\textsc{FP}\times\textsc{FN}}{\sqrt{(\textsc{TN}+\textsc{FN})(\textsc{TN}+\textsc{FP})(\textsc{TP}+\textsc{FP})(\textsc{TP}+\textsc{FN})}}$

where TP is the number of true positives, TN the number of true negatives, FP
the number of false positives and FN the number of false negatives. We also
provide the precision and recall scores for they are more intuitive. Precision
is a proportion of relevant instances among retrieved instances, while recall
is the proportion of relevant instances that have been retrieved. They are
defined as:

$\textsc{Precision}=\frac{\textsc{TP}}{\textsc{TP}+\textsc{FP}}\\\
\textsc{Recall}=\frac{\textsc{TP}}{\textsc{TP}+\textsc{FN}}$

The classical decision rules and the result of their application to the test
dataset after augmentation and imputation are summarized in Tab. 2. The MCC
score was calculated on the test part of the dataset. It should be mentioned,
that even though the filters in different surveys have similar names (i.e.
$y_{PS1}$ and $Y_{DES}$), they are not identical to each other. Therefore, it
is not entirely correct to apply the rules made for one survey to the
magnitudes of the other survey. However, we have estimated that for our
dataset the magnitudes of the same name differ within $0.2$ mag which does not
increase the score in any case. Also, it is worth mentioning that Burningham
et al. (2013) was originally devoted to T-type dwarfs exclusively, yet it
shows great performance on L dwarfs, so we use it as a decision rule for both
L and T dwarfs.

Although the performance of the decision rules is reasonably high, the actual
number of false positive and false negative classifications grows with the
number of objects, and this becomes important when we have millions of
objects, as in most modern sky surveys (PanSTARRS - 1.9 billion objects, 2MASS
- 470 million objects, ALLWISE - 560 million objects), so it makes sense to
make an effort to increase the performance and these values (Tab. 2) are the
baselines we want to outperform.

Table 2: Decision rules from the literature. Author | Rule | MCC
---|---|---
Carnero Rosell et al. (2019) |  $(i-z)>1.2$, $(z-Y)>0.15$, | 0.935
| ($Y_{AB}-J_{Vega})>1.6$, z¡ 22 |
Burningham et al. (2013) |  $(z-J)_{Vega}>2.5$, $J<17.5$ | 0.921

Figure 5: The confusion matrices for the decision rules Carnero Rosell et al.
(2019) on the top and Burningham et al. (2013) on the bottom applied to the
test part of the dataset.

The work made use of the following Python software packages: Sckit-learn 1.0.2
version, optuna 2.10.1 version, Tabnet from Pytorch 4.0 version.

###  3.1 Random Forest

The decision tree concept is quite similar to the classical colour selections
that are traditionally used in astronomy to classify objects. While automated
decision trees can be much more efficient than classic decision rules, they
tend to overfit, i.e. they learn too much about the data they are trained on
and can fail when applied to data they have not seen before. The solution to
this problem can be a random forest approach (RF) - an ensemble of decision
trees. In this case, the decision about the class an object belongs to, is
made based on which class the greater number of trees voted for.

Using optuna, we selected the maximum tree depth, the minimum number of
samples required to split an internal node, the criterion and the maximum
number of features in a node. Tuned parameters and corresponding scores are
presented in Tab. 3. Note that the scores are not guaranteed and depend on the
dataset the model is tested on, e.g. the imputing of missing values. For a
more objective evaluation of the model we use bootstrap technique (see
confidence intervals, Sec. 4 and Fig. 9).

In Fig. 6, an example of a two-dimensional cut shows the separating boundary
between classes defined by the RF model. It should be noted that this is a
slice in multidimensional parameter space, so it does not represent the number
of mismatches. In fact, the performance is very good, meaning that the
decision rules underlying that model might be quite complex and cannot be
represented on a 2D diagram.

Fig. 10 shows the importance of features for the RF model in all three cases.
Seeing Pan-STARRS magnitudes, the model primarily relies on $i$ magnitudes.
Without optical magnitudes, $i-y$ becomes the most important feature. This
becomes even more pronounced for the RF model that relies on colours only.

Figure 6: Slice of the separating boundary in the feature space according to the RF model. Table 3: Random forest hyperparameters for three sets of features. Number of trees is 500 for all of the models. Number in max features is the fraction of all available features. Hyperparameter | All features | No PS magnitudes | Only colours
---|---|---|---
Max depth | 11 | 13 | 12
Min samples split | 20 | 9 | 8
Max features | sqrt | sqrt | 0.7
Criterion | entropy | entropy | gini
MCC test score | 0.983 | 0.986 | 0.975
MCC train score | 0.987 | 0.990 | 0.990
Precision | 0.992 | 0.992 | 0.988
Recall | 0.992 | 0.994 | 0.987

###  3.2 XGBoost

Boosting is a very popular machine-learning technique. It is a type of
ensemble learning that uses the output of the previous model as input to the
next one. Instead of training models individually, boosting trains models
sequentially, with each new model trained to correct the errors of previous
ones. At each iteration, correctly predicted results are given less weight,
and incorrectly predicted results are given more weight. It then uses the
weighted average to get the final result.

We are also interested in boosting from the point of view that often a certain
principle of handling missing values is built into the algorithm. Two popular
boosting models are CatBoost (Prokhorenkova et al., 2017) and XGBoost (Chen
and Guestrin, 2016). CatBoost only has built-in filling in the missing values
with some specific values, but XGBoost uses the following strategy: each node
is assigned a default solution, and this works well in many cases. Therefore,
we choose XGBoost for the task and also compare the performance of the model
on the missing values filled in by the default method of XGBoost and filling
with the Iterative Imputer.

On a test dataset with the default missing value algorithm, XGBoost gives
$\textsc{MCC}=0.96$. When training and testing on data in which we filled the
missing values using the Iterative Imputer, the performance reaches
$\textsc{MCC}=0.986$. Thus, we conclude that Iterative Imputer in this case is
not only a more robust method but also has a better effect on model
performance.

The number of estimators is fixed at 500. The following hyperparameters were
optimized with optuna: maximum tree depth, learning rate (the step size
reduction used during the update to prevent overfitting), subsample ratio of
the training instances, reg_alpha (regularization term) and gamma - the
minimum loss reduction required to create a further partition on a leaf node
of the tree. Optimized values are presented in Tab. 4.

Table 4: XGBoost hyperparameters for three sets of features. Hyperparameter | All features | No PS magnitudes | Only colours
---|---|---|---
Max depth | 15 | 15 | 10
Learning rate | 0.340 | 0.126 | 0.033
Subsample | 0.05 | 0.04 | 0.11
Gamma | 0.323 | 0.548 | 0.996
Reg alpha | 0.82 | 0.48 | 0.03
MCC test score | 0.978 | 0.972 | 0.969
MCC train score | 0.980 | 0.978 | 0.974
Precision | 0.992 | 0.989 | 0.985
Recall | 0.985 | 0.982 | 0.985

The feature importances is presented in Fig. 10. Trained on all features
XGBoost relies mainly on the $i$ magnitude of the Pan-STARRS survey, as does
RF. If Pan-STARRS magnitudes are excluded from the feature list, $i-z$, $z-J$
and $i-y$ colours became the dominant features, however, the $J$ magnitude of
the 2MASS survey remains important.

###  3.3 Support Vector Machine

The support vector machine (Cortes and Vapnik, 1995) is another widely used
and well-developed method. The principle of the support vector machine is to
find a line, surface or hypersurface that would separate classes in the
feature space. The fitting process maximizes the distance from each point to
the decision boundary (reference vector).

Figure 7: Slice of the separating boundary in the space of features according
to the SVM model.

We tuned the regularization parameter C, the kernel type and the kernel
coefficient (’gamma’, for ’rbf’ kernel). The decision function was set to one-
versus-one (’ovo’) since it is a binary classification and weights of classes
were automatically adjusted inversely proportional to class frequencies in the
input data. For tuned parameters in every case see Tab. 5. It is seen from the
table that models are very similar for all of the cases and so are the most
important features (see Fig. 10). As SVM mainly rely on colour indices, the
importance distributions do not change significantly when some or all of the
magnitudes get excluded.

Table 5: SVM classifier hyperparameters for three sets of features. Hyperparameter | All features | No PS magnitudes | Only colours
---|---|---|---
Kernel | rbf | linear | rbf
C | 1.150 | 0.729 | 0.792
Gamma | 0.298 | 0.066 | 0.453
MCC test score | 0.981 | 0.984 | 0.958
MCC train score | 0.982 | 0.980 | 0.968
Precision | 0.989 | 0.986 | 0.972
Recall | 0.992 | 0.998 | 0.986

Fig. 7 shows a separating boundry constructed by the SVM model using a 2D cut
as an example. It is worth noting that this is only a slice in which the
remaining feature values are taken in some neighbourhood of the mean, so a
large number of points that fell into the wrong area does not mean an actual
misclassification.

###  3.4 TabNet

TabNet (Arik and Pfister, 2019) is a deep learning neural network that
utilises attention to select important features at each step of the decision-
making process so that only the most important features are used. In this
case, the choice of features depends on the object, and, for example, it can
be different for each row of the training data set. In the end, you can see
which features the model has focused on the most.

TabNet consists of several steps, each step is a block of components, with the
number of steps being a hyperparameter. Each step gets its vote in the final
classification, which mimics the ensemble classification. Other
hyperparameters are the width of the decision prediction layer (N_d), the
width of the attention embedding for each mask (N_a), number of shared Gated
Linear Units at each step (N_shared) and the coefficient for feature reusage
in the masks (Gamma). We fitted the hyperparameters of TabNet on optuna using
MCC score as the metric. The fitted parameters are presented in Tab. 6.

Table 6: TabNet hyperparameters for the different sets of features. Hyperparameter | All features | No PS magnitudes | Only colours
---|---|---|---
N_d | 16 | 12 | 12
N_a | 16 | 28 | 12
Number of steps | 3 | 3 | 2
Gamma | 1.2 | 1.6 | 1.6
N_shared | 2 | 2 | 1
MCC test score | 0.983 | 0.986 | 0.975
MCC train score | 0.987 | 0.990 | 0.990
Precision | 0.992 | 0.992 | 0.988
Recall | 0.992 | 0.994 | 0.987

The model is trained using gradient descent optimization with Adam as the
optimizer. The validation part of a dataset is used in training to prevent
overfitting, so the result of the model training is the configuration, that
provides the best scores, both on training and validation data. A separating
boundary is shown in Fig. 8, as it can be seen it is more complex then the
boundaries of other models.

Figure 8: Slice of the separating boundary in the space of features according
to the TabNet model.

##  4 Results and discussion

In this section, we summarise the results of the application of machine
learning methods to the brown dwarf search problem.

We trained four models on the data: Random Forest Classifier, SVM Classifier,
XGBoost Classifier and TabNet Classifier.

Confidence intervals for the MCC scores of the obtained models are presented
in Fig. 9. Confidence intervals are calculated via the bootstrap method with
100 samples half the length of the test data set. The coloured box is the
interval from the 25th percentile to the 75th percentile and the median value
is represented by a black line. Error bars show the minimum and the maximum
value with the outliers marked as diamonds. The baseline values obtained using
decision rules from the literature are represented in Fig. 9 as dashed lines.

All models on the full set of features and the features without PS1 magnitudes
provide almost the same results, but some performed slightly better. If we use
only magnitude features, the performance decreases. However, such models still
outperform the baselines.

Although the performance of the models is nearly the same, they differ in
terms of robustness. The importance of features of the models is presented in
Fig. 10. Random Forest and XGBoost rely primarily on the $i$ magnitude in
their decisions, if any. At the same time, SVM and Tabnet on the full set of
features seem more robust, as they rely mainly on colour indices. Tabnet also
tends to pay a lot of attention to the magnitudes, although it performs the
best in the case of using only colours. Table 7 shows true positive, true
negative, false positive and false negative values for all of the models
obtained on the test part of the dataset.

According to Carnero Rosell et al. (2019) and Skrzypek et al. (2015),
$(i-z)_{PS1}$ and $y_{PS1}-J$ colour indices are expected to be the most
important feature since they have the largest variation across the M/L
transition. While the $(i-z)_{PS1}$ colour index is important to SVM
classifiers and in some cases for RF and TabNet, most of the models do not
consider them essential.

It is expected that $z_{PS1}-J$ is the most important feature as Burningham et
al. (2013), which we used as a baseline, is nearly entirely based on the above
colour. However, it only plays a secondary role in most of the models. Unlike
the previous works, we revealed the importance of the $(i-y)_{PS1}$ colour
index. It is the most important feature in most cases and the colour selection
$(i-y)_{PS1}>1.88$ alone gives an MCC score of $0.968$ on the testing data.
Other colour indices could be important in the case of multi-class
classification, for example, L and T-type dwarfs differ drastically in
$W1-W2$, although it has almost no relevance for the problem investigated in
the present article.

Table 7: True positive (TP), true negative (TN), false positive (FP) and false negative (FN) values as well as Precison and Recall scores for four models: Random Forest (RF), XGBoost, Support Vector Machine (SVM), and TabNet. Each model was trained on three sets of features, labeled as ”All features”, ”w/o PS magnitudes”, and ”only colours”. Model | TP | TN | FP | FN | Precision | Recall
---|---|---|---|---|---|---
Random Forest
All features | 846 | 7 | 846 | 7 | 0.992 | 0.992
W/o PS magnitudes | 848 | 7 | 833 | 5 | 0.992 | 0.994
Only colours | 842 | 10 | 830 | 11 | 0.988 | 0.987
XGBoost
All features | 840 | 6 | 834 | 13 | 0.992 | 0.985
W/o PS magnitudes | 838 | 9 | 831 | 15 | 0.989 | 0.982
Only colours | 840 | 13 | 827 | 13 | 0.985 | 0.985
SVM
All features | 846 | 9 | 831 | 7 | 0.989 | 0.992
W/o PS magnitudes | 851 | 12 | 828 | 2 | 0.986 | 0.998
Only colours | 841 | 24 | 816 | 12 | 0.972 | 0.986
TabNet
All features | 850 | 9 | 831 | 3 | 0.992 | 0.992
W/o PS magnitudes | 851 | 12 | 828 | 2 | 0.992 | 0.994
Only colours | 846 | 13 | 827 | 7 | 0.988 | 0.987
Figure 9: Confidence intervals for the scores of models Figure 10: Importance
of features for all models. For RF, XGBoost and SVM models we calculate the
importance of each feature using SHAP and for TabNet we use built-in
capabilities of the model.

##  5 Conclusion

In this paper, we compiled a dataset of L and T type brown dwarfs (marked as
positive class) and objects of other spectral types (marked as negative class)
from literature. It is important to acknowledge potential limitations and
biases that could impact the dataset’s representativeness and
generalizability. While we have compiled the data set in such a way as to
reproduce the observed distribution in absolute magnitudes, the distribution
in apparent magnitudes may be not representative, shifted towards brighter
objects. This most likely does not affect the performance of the model in the
case of ”only colours” mode, but may have an impact on models trained on
magnitudes too. The incompleteness in M type dwarfs that stems from
limitations in the Best et al. (2018) catalogue introduces a bias in the
dataset towards the included objects, potentially impacting the
generalizability of the results.

Brown dwarfs are faint astronomical sources with peak intensities that fall
into the infra-red part of the spectrum, so it is harder to carry out the
measurements in optical and far-red filters, such as $i$ and $z$. We imputed
missing values with Iterative Imputer and explored the result. We also imputed
the colour indices regardless of corresponding magnitudes in order to reduce
the error of imputing. For most of the magnitude features, the imputing error
is compatible with the measurement error. For colour index features, imputing
error is usually much lower than the corresponding measurement error. Thus,
the imputing part of preprocessing is considered to be successful.

Four models, namely, Random Forest Classifier, SVM Classifier, XGBoost
Classifier and TabNet Classifier were trained to distinguish brown dwarfs
among all objects according to their photometric data. The classification
results for all models are consistently high, all models outperform the
baselines. However, tree-like models (RF and XGBoost) tend to exploit the
faintness of the objects, which is less preferable in terms of robustness of
the model. On the contrary, SVM primarily relies on colour indices, which
reduce the possibility of misclassification of other types of faint sources.

Examining the features of the models, we have found that the $(i-y)_{PS1}$
colour index is the most important feature in most cases. This colour index
can be used subsequently and independently in brown dwarf classification
problems. We also confirm that the $z_{PS1}-J$, $(i-z)_{PS1}$ and
$(z-y)_{PS1}$ colour indices are important and powerful features in brown
dwarf classification.

Along with the overall success of the work, there are several limitations to
the applicability of the results, that are yet to be solved. First, potential
biases of the dataset were emphasised, including a bias towards brighter
objects. These biases can affect the generalizability and reliability of the
findings. Second, it is known, that highly red-shifted quasars might also be a
source of contamination when distinguishing brown dwarfs from other types of
objects (Carnero Rosell et al., 2019). This issue can be resolved using colour
cuts in most cases, but can still remain for some cases (Reed et al., 2017).
At this point, we haven’t considered quasars. Third, while the imputing
magnitudes and colour indices independently allow us to predict colour indices
much more accurately, this leads to the colour index not being always a linear
combination of magnitudes involved. Multi-class classification will also be in
scope of future work.

## Acknowledgements

This work was supported by by Non-commercial Foundation for the Advancement of
Science and Education INTELLECT.

Author is grateful to Ilya Dyugay, Konstantin Malanchev, Dana Kovaleva and
Oleg Malkov for the fruitful discussion and valuable advice. Author thanks the
anonymous referee for careful reading and suggestions, which greatly helped to
improve the article.

This research has made use of NASA’s Astrophysics Data System Bibliographic
Services.

## References

  * Kumar (1963) Kumar, S.S.,  The structure of stars of very low mass.,  The Astrophysical Journal 137 (1963) 1121.
  * Hayashi and Nakano (1963) Hayashi, C. and Nakano, T.,  Evolution of stars of small masses in the pre-main-sequence stages,  Progress of Theoretical Physics 30 (1963) 460–474.
  * Rebolo et al. (1995) Rebolo, R., Zapatero Osorio, M.R., and Martín, E.L.,  Discovery of a brown dwarf in the pleiades star cluster,  Nature 377 (1995) 129–131.
  * Nakajima et al. (1995) Nakajima, T., Oppenheimer, B.R., Kulkarni, S.R., et al.,  Discovery of a cool brown dwarf,  Nature 378 (1995) 463–465.
  * Luhman (2013) Luhman, K.L.,  Discovery of a binary brown dwarf at 2 pc from the sun,  The Astrophysical Journal Letter 767 (2013) L1.
  * Burningham et al. (2013) Burningham, B., Cardoso, C.V., Smith, L., et al.,  76 t dwarfs from the ukidss las: benchmarks, kinematics and an updated space density,  Monthly Notices of the Royal Astronomical Society 433 (2013) 457–497.
  * Carnero Rosell et al. (2019) Carnero Rosell, A., Santiago, B., dal Ponte, M., et al.,  Brown dwarf census with the dark energy survey year 3 data and the thin disc scale height of early l types,  Monthly Notices of the Royal Astronomical Society 489 (2019) 5301–5325.
  * Kirkpatrick et al. (1999) Kirkpatrick, J.D., Reid, I.N., Liebert, J., et al.,  Dwarfs Cooler than “M“: The Definition of Spectral Type “L” Using Discoveries from the 2 Micron All-Sky Survey (2MASS),  The Astrophysical Journal 519 (1999) 802–833.
  * Skrzypek et al. (2016) Skrzypek, N., Warren, S.J., and Faherty, J.K.,  Vizier online data catalog: Photometric brown-dwarf classification (skrzypek+, 2016),  VizieR Online Data Catalog (2016) J/A+A/589/A49.
  * Kirkpatrick et al. (2021) Kirkpatrick, J.D., Gelino, C.R., Faherty, J.K., et al.,  The field substellar mass function based on the full-sky 20 pc census of 525 l, t, and y dwarfs,  Astrophysical Journal Supplement Series 253 (2021) 7.
  * Mužić et al. (2017) Mužić, K., Schödel, R., Scholz, A., et al.,  The low-mass content of the massive young star cluster rcw 38,  Monthly Notices of the Royal Astronomical Society 471 (2017) 3699–3712.
  * Smith et al. (2014) Smith, L., Lucas, P.W., Bunce, R., et al.,  High proper motion objects from the UKIDSS Galactic plane survey,  Monthly Notices of the Royal Astronomical Society 443 (2014) 2327–2341.
  * Lodieu et al. (2014) Lodieu, N., Pérez-Garrido, A., Béjar, V.J.S., et al.,  Binary frequency of planet-host stars at wide separations. A new brown dwarf companion to a planet-host star,  Astronomy and Astrophysics 569 (2014) A120.
  * Artigau et al. (2009) Artigau, É., Bouchard, S., Doyon, R., et al.,  Photometric variability of the t2.5 brown dwarf simp j013656.5+093347: Evidence for evolving weather patterns,  The Astrophysical Journal 701 (2009) 1534 – 1539.
  * Gillon et al. (2013) Gillon, M., Triaud, A.H.M.J., Jehin, E., et al.,  Fast-evolving weather for the coolest of our two new substellar neighbours,  Astronomy & Astrophysics 555 (2013) L5.
  * Khandrika et al. (2013) Khandrika, H.G., Burgasser, A.J., Melis, C., et al.,  A search for photometric variability in l- and t-type brown dwarf atmospheres,  The Astronomical Journal 145 (2013).
  * Marley and Robinson (2015) Marley, M.S. and Robinson, T.D.,  On the cool side: Modeling the atmospheres of brown dwarfs and giant planets,  Annual Review of Astronomy and Astrophysics 53 (2015) 279–323.
  * Charnay et al. (2018) Charnay, B., Bézard, B., Baudino, J.L., et al.,  A self-consistent cloud model for brown dwarfs and young giant exoplanets: Comparison with photometric and spectroscopic observations,  The Astrophysical Journal 854 (2018) 172.
  * Tremblin et al. (2019) Tremblin, P., Padioleau, T., Phillips, M.W., et al.,  Thermo-compositional diabatic convection in the atmospheres of brown dwarfs and in earth’s atmosphere and oceans,  The Astrophysical Journal 876 (2019) 144.
  * Tan and Showman (2019) Tan, X. and Showman, A.P.,  Atmospheric variability driven by radiative cloud feedback in brown dwarfs and directly imaged extrasolar giant planets,  The Astrophysical Journal 874 (2019) 111.
  * Burningham et al. (2017) Burningham, B., Marley, M.S., Line, M.R., et al.,  Retrieval of atmospheric properties of cloudy l dwarfs,  Monthly Notices of the Royal Astronomical Society 470 (2017) 1177–1197.
  * Saumon and Marley (2008) Saumon, D. and Marley, M.S.,  The evolution of l and t dwarfs in color-magnitude diagrams,  The Astrophysical Journal 689 (2008) 1327–1344.
  * Vos et al. (2019) Vos, J.M., Allers, K., Apai, D., et al.,  Astro2020 white paper: The l/t transition (2019).
  * Maravelias et al. (2022) Maravelias, G., Bonanos, A.Z., Tramper, F., et al.,  A machine-learning photometric classifier for massive stars in nearby galaxies i. the method,  arXiv e-prints (2022) arXiv:2203.08125.
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al.,  Scikit-learn: Machine learning in Python,  Journal of Machine Learning Research 12 (2011) 2825–2830.
  * Mohammadi et al. (2022) Mohammadi, M., Mutatiina, J., Saifollahi, T., et al.,  Detection of extragalactic ultra-compact dwarfs and globular clusters using explainable AI techniques,  Astronomy and Computing 39 (2022) 100555.
  * Best et al. (2018) Best, W.M.J., Magnier, E.A., Liu, M.C., et al.,  Photometry and proper motions of m, l, and t dwarfs from the pan-starrs1 3$\pi$ survey,  Astrophysical Journal Supplement Series 234 (2018) 1.
  * Chambers et al. (2016) Chambers, K.C., Magnier, E.A., Metcalfe, N., et al.,  The pan-starrs1 surveys,  arXiv e-prints (2016) arXiv:1612.05560.
  * Cutri et al. (2003) Cutri, R.M., Skrutskie, M.F., van Dyk, S., et al.,  Vizier online data catalog: 2mass all-sky catalog of point sources (cutri+ 2003),  VizieR Online Data Catalog (2003) II/246.
  * Cutri et al. (2021) Cutri, R.M., Wright, E.L., Conrow, T., et al.,  Vizier online data catalog: Allwise data release (cutri+ 2013),  VizieR Online Data Catalog (2021) II/328.
  * Collaboration et al. (2016) Collaboration, G., Prusti, T., de Bruijne, J.H.J., et al.,  The gaia mission,  Astronomy and Astrophysics 595 (2016) A1.
  * Collaboration et al. (2022) Collaboration, G., Vallenari, A., Brown, A.G.A., et al.,  Gaia data release 3: Summary of the content and survey properties,  arXiv e-prints (2022) arXiv:2208.00211.
  * Akiba et al. (2019) Akiba, T., Sano, S., Yanase, T., et al.,  Optuna: A next-generation hyperparameter optimization framework (2019).
  * Lundberg and Lee (2017) Lundberg, S.M. and Lee, S.I.,  A unified approach to interpreting model predictions,  in: Guyon, I., Luxburg, U.V., Bengio, S., et al. (Eds.), Advances in Neural Information Processing Systems 30, Curran Associates, Inc., 2017, pp. 4765–4774. URL: http://papers.nips.cc/paper/7062-a-unified-approach-to-interpreting-model-predictions.pdf.
  * Prokhorenkova et al. (2017) Prokhorenkova, L., Gusev, G., Vorobev, A., et al.,  CatBoost: unbiased boosting with categorical features,  arXiv e-prints (2017) arXiv:1706.09516.
  * Chen and Guestrin (2016) Chen, T. and Guestrin, C.,  Xgboost: A scalable tree boosting system (2016).
  * Cortes and Vapnik (1995) Cortes, C. and Vapnik, V.,  Support-vector networks,  Machine learning 20 (1995) 273–297.
  * Arik and Pfister (2019) Arik, S.O. and Pfister, T.,  Tabnet: Attentive interpretable tabular learning,  arXiv e-prints (2019) arXiv:1908.07442.
  * Skrzypek et al. (2015) Skrzypek, N., Warren, S.J., Faherty, J.K., et al.,  Photometric brown-dwarf classification,  Astronomy & Astrophysics 574 (2015) A78.
  * Reed et al. (2017) Reed, S.L., McMahon, R.G., Martini, P., et al.,  Eight new luminous z $\geq$ 6 quasars discovered via sed model fitting of vista, wise and dark energy survey year 1 observations,  Monthly Notices of the Royal Astronomical Society 468 (2017) 4702–4718.


# Stability Analysis of Higher-order Ecological Networks
through Tensor Decompositions

Anqi Dong1 and Can Chen2 1Anqi Dong is with the Department of Mechanical and
Aerospace Engineering, University of California, Irvine, Irvine, CA 92617,
USA<EMAIL_ADDRESS>Chen is with the School of Data Science and Society
and the Department of Mathematics, University of North Carolina at Chapel
Hill, Chapel Hill, NC 27599, USA<EMAIL_ADDRESS>
###### Abstract

Complex ecological networks are characterized by intricate interactions that
extend beyond traditional pairwise relationships. Understanding the stability
of higher-order ecological networks is salient for species coexistence,
biodiversity, and community persistence. In this article, we present a novel
approach for stability analysis of higher-order ecological networks through
tensor decompositions. We are particularly 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.

###### Index Terms:

Linear stability, Jacobian matrices, ecological networks, higher-order
interactions, tensor decompositions.

##  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 [31, 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) 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 [23, 15,
28], have found applications across a wide array of fields, including
dynamical systems [11, 13, 9], signal processing [29, 14], network science
[10, 12, 32], data analysis [33, 24], and more.

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) [26, 2, 19]. 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 [21, 30, 20].

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) [23] 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. This is because CPD is not
numerically stable. Moreover, the low-rank approximation of CPD is ill-posed.

  * •

The tensor train decomposition (TTD) [28] 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. More 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 is
often characterized by the higher-order generalized Lotka–Volterra (HOGLV)
model [5, 1, 30, 20], 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 given by

$\begin{split}\dot{\textbf{x}}&=\textbf{x}*\Big{[}\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\Big{]}.\end{split}$

Remark 1: If the maximum order of interactions is $M$, the total number of the
model parameters in the HOGLV model can be estimated as

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

Obviously, the above 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
follows:

$\begin{split}\mathcal{B}&\approx\mathcal{B}_{0}\times_{123}\\{\textbf{B}_{1},\textbf{B}_{2},\textbf{B}_{3}\\},\\\
\mathcal{C}&\approx\mathcal{C}_{0}\times_{1234}\\{\textbf{C}_{1},\textbf{C}_{2},\textbf{C}_{3},\textbf{C}_{4}\\},\end{split}$

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 2: 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 also 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

Second, 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
follows:

$\begin{split}\mathcal{B}&\approx\sum_{i=1}^{b}[\textbf{B}_{1}]_{:i}\circ[\textbf{B}_{2}]_{:i}\circ[\textbf{B}_{3}]_{:i},\\\
\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},\\\
\end{split}$

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 3: 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

Lastly, 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
follows:

$\begin{split}\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}},\\\
\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}},\end{split}$

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 4: 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. Most importantly, TTD can offer the crucial
advantages of numerical stability and quasi-optimal approximation.

##  IV 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, 26]. 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. In
this section, 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. ∎

Remark 5: 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 6: 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 7: 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 8: 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 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 [25, 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 [22] and TT Toolbox [27].

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

First, 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 defined as
$\textbf{x}^{*}=\begin{bmatrix}1&1&1&1\end{bmatrix}^{\top}$ for both networks,
the associated Jacobian matrices can be computed as

$\displaystyle\textbf{M}_{1}$
$\displaystyle=\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},$
$\displaystyle\textbf{M}_{2}$
$\displaystyle=\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

Next, 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

Lastly, 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 of the HOGLV model.

##  VI Conclusion

To summarize, 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 matrices and tenors are time-varying, is a promising direction for
future research.

## References

  * [1] Mohammad AlAdwani and Serguei Saavedra.  Is the addition of higher-order interactions in ecological models increasing the understanding of ecological dynamics?  Mathematical Biosciences, 315:108222, 2019.
  * [2] S. Allesina and S. Tang.  Stability criteria for complex ecosystems.  Nature, 483:205–208, 2012.
  * [3] Stefano Allesina and Si Tang.  Stability criteria for complex ecosystems.  Nature, 483(7388):205–208, 2012.
  * [4] Stefano Allesina and Si Tang.  The stability–complexity relationship at age 40: a random matrix perspective.  Population Ecology, 57(1):63–75, 2015.
  * [5] Eyal Bairey, Eric D. Kelsic, and Roy Kishony.  High-order species interactions shape ecosystem diversity.  Nature communications, 7(1):12285, 2016.
  * [6] Joseph W. Baron, Thomas Jun Jewell, Christopher Ryder, and Tobias Galla.  Breakdown of random-matrix universality in persistent lotka-volterra communities.  Physical Review Letters, 130(13):137401, 2023.
  * [7] Subhash C. Bhargava.  Generalized lotka-volterra equations and the mechanism of technological substitution.  Technological Forecasting and Social Change, 35(4):319–326, 1989.
  * [8] Guy Bunin.  Ecological communities with lotka-volterra dynamics.  Physical Review E, 95(4):042414, 2017.
  * [9] Can Chen.  Explicit solutions and stability properties of homogeneous polynomial dynamical systems.  IEEE Transactions on Automatic Control, 68(8):4962–4969, 2023.
  * [10] Can Chen and Indika Rajapakse.  Tensor entropy for uniform hypergraphs.  IEEE Transactions on Network Science and Engineering, 7(4):2889–2900, 2020.
  * [11] Can Chen, Amit Surana, Anthony Bloch, and Indika Rajapakse.  Multilinear time invariant system theory.  In 2019 Proceedings of the Conference on Control and its Applications. SIAM, 2019.
  * [12] Can Chen, Amit Surana, Anthony M. Bloch, and Indika Rajapakse.  Controllability of hypergraphs.  IEEE Transactions on Network Science and Engineering, 8(2):1646–1657, 2021.
  * [13] Can Chen, Amit Surana, Anthony M. Bloch, and Indika Rajapakse.  Multilinear control systems theory.  SIAM Journal on Control and Optimization, 59(1):749–776, 2021.
  * [14] Andrzej Cichocki, Danilo Mandic, Lieven De Lathauwer, Guoxu Zhou, Qibin Zhao, Cesar Caiafa, and Huy Anh Phan.  Tensor decompositions for signal processing applications: From two-way to multiway component analysis.  IEEE signal processing magazine, 32(2):145–163, 2015.
  * [15] Lieven De Lathauwer, Bart De Moor, and Joos Vandewalle.  A multilinear singular value decomposition.  SIAM journal on Matrix Analysis and Applications, 21(4):1253–1278, 2000.
  * [16] Sergey V. Dolgov, Boris N. Khoromskij, Ivan V. Oseledets, and Dmitry V. Savostyanov.  Computation of extreme eigenvalues in higher dimensions using block tensor train format.  Computer Physics Communications, 185(4):1207–1216, 2014.
  * [17] Anqi Dong, Tryphon T. Georgiou, and Allen Tannenbaum.  Data assimilation for sign-indefinite priors: A generalization of sinkhorn’s algorithm.  arXiv preprint arXiv:2308.11791, 2023.
  * [18] Anqi Dong, Tryphon T. Georgiou, and Allen Tannenbaum.  Negative probabilities in gene regulatory networks.  arXiv preprint arXiv:2307.07738, 2023.
  * [19] Theo Gibbs, Jacopo Grilli, Tim Rogers, and Stefano Allesina.  Effect of population abundances on the stability of large random ecosystems.  Physical Review E, 98(2):022410, 2018.
  * [20] Theo Gibbs, Simon A. Levin, and Jonathan M. Levine.  Coexistence in diverse communities with higher-order interactions.  Proceedings of the National Academy of Sciences, 119(43):e2205063119, 2022.
  * [21] Jacopo Grilli, György Barabás, Matthew J Michalska-Smith, and Stefano Allesina.  Higher-order interactions stabilize dynamics in competitive network models.  Nature, 548(7666):210–213, 2017.
  * [22] Tamara Kola, Brett W. Bader, Evrim NMN Acar Ataman, Daniel Dunlavy, Robert Bassett, Casey J. Battaglino, Todd Plantenga, Eric Chi, and Samantha Hansen.  Tensor toolbox for matlab.  Technical report, Sandia National Lab, Albuquerque, NM (United States), 2017.
  * [23] Tamara G. Kolda and Brett W. Bader.  Tensor decompositions and applications.  SIAM review, 51(3):455–500, 2009.
  * [24] Yipeng Liu, Jiani Liu, Zhen Long, and Ce Zhu.  Tensor computation for data analysis.  Springer, 2022.
  * [25] Thomas Mach.  Computing inner eigenvalues of matrices in tensor train matrix format.  In Numerical Mathematics and Advanced Applications 2011: Proceedings of the 9th European Conference on Numerical Mathematics and Advanced Applications, Leicester, September 2011, pages 781–788. Springer, 2012.
  * [26] Robert M May.  Will a large complex system be stable?  Nature, 238(5364):413–414, 1972.
  * [27] IV Oseledets, S. Dolgov, V. Kazeev, D. Savostyanov, O Lebedeva, P Zhlobich, T. Mach, and L. Song.  TT-toolbox.  URL: http://oseledets. github. io/software, 2016.
  * [28] Ivan V. Oseledets.  Tensor-train decomposition.  SIAM Journal on Scientific Computing, 33(5):2295–2317, 2011.
  * [29] Nicholas D Sidiropoulos, Lieven De Lathauwer, Xiao Fu, Kejun Huang, Evangelos E Papalexakis, and Christos Faloutsos.  Tensor decomposition for signal processing and machine learning.  IEEE Transactions on signal processing, 65(13):3551–3582, 2017.
  * [30] Pragya Singh and Gaurav Baruah.  Higher order interactions and species coexistence.  Theoretical Ecology, 14:71–83, 2021.
  * [31] Per Sebastian Skardal and Alex Arenas.  Higher order interactions in complex networks of phase oscillators promote abrupt synchronization switching.  Communications Physics, 3(1):218, 2020.
  * [32] Amit Surana, Can Chen, and Indika Rajapakse.  Hypergraph similarity measures.  IEEE Transactions on Network Science and Engineering, 10(2):658–674, 2022.
  * [33] Hua Zhou, Lexin Li, and Hongtu Zhu.  Tensor regression with applications in neuroimaging data analysis.  Journal of the American Statistical Association, 108(502):540–552, 2013.