# MulGuisin, a Topological Network Finder and its Performance on Galaxy
Clustering

Young Ju1,2, Inkyu Park1,2, Cristiano G. Sabiu1,2 and Sungwook E. Hong3,4
1Department of Physics, University of Seoul, 163 Seoulsiripdae-ro, Dongdaemun-
gu, Seoul 02504, Republic of Korea
2Natural Science Research Institute, University of Seoul, 163 Seoulsiripdae-
ro, Dongdaemun-gu, Seoul 02504, Republic of Korea
3Korea Astronomy and Space Science Institute, 776 Daedeok-daero, Yuseong-gu,
Daejeon 34055, Republic of Korea
4Astronomy Campus, University of Science and Technology, 776 Daedeok-daero,
Yuseong-gu, Daejeon 34055, Republic of Korea E-mail<EMAIL_ADDRESS>
###### Abstract

We introduce a new clustering algorithm, MulGuisin (MGS), that can identify
distinct galaxy over-densities using topological information from the galaxy
distribution. This algorithm was first introduced in an LHC experiment as a
Jet Finder software, which looks for particles that clump together in close
proximity. The algorithm preferentially considers particles with high energies
and merges them only when they are closer than a certain distance to create a
jet. MGS shares some similarities with the minimum spanning tree (MST) since
it provides both clustering and network-based topology information. Also,
similar to the density-based spatial clustering of applications with noise
(DBSCAN), MGS uses the ranking or the local density of each particle to
construct clustering. In this paper, we compare the performances of clustering
algorithms using controlled data and some realistic simulation data as well as
the SDSS observation data, and we demonstrate that our new algorithm find
networks most efficiently and it defines galaxy networks in a way that most
closely resembles human vision.

###### keywords:

large-scale structure of Universe, galaxies: clusters: general, methods:
statistical, software: data analysis

††pagerange: MulGuisin, a Topological Network Finder and its Performance on
Galaxy Clustering–MulGuisin, a Topological Network Finder and its Performance
on Galaxy Clustering

##  1 Introduction

In the standard $\Lambda$CDM cosmology paradigm, structures in the universe
grow in a hierarchical manner (e.g., White & Rees, 1978; Fall & Efstathiou,
1980; Blumenthal et al., 1984). It means that smaller structures of matter
start forming earlier, and the more massive structures form by the merging and
accretion of smaller structures at later epoch of the universe. Therefore,
understanding cosmic structures in the universe at various scales is crucial
for understanding the nature of our universe. For example, numerous statistics
of large-scale structures, such as topological analyses (e.g., Gott et al.,
1986; Park & Gott, 1991; Park & Kim, 2010; Appleby et al., 2017; Appleby et
al., 2018), Alcock-Paczynski tests (e.g., Alcock & Paczynski, 1979; Ballinger
et al., 1996; Li et al., 2014; Park et al., 2019), and small-scale redshift-
space distortion (RSD; e.g., Sheth, 1996; DeRose et al., 2019; Tonegawa et
al., 2020), have been used to constraint cosmological parameters such as the
matter density parameter ($\Omega_{\rm m}$) and the equation-of-state
parameter of dark energy ($\omega_{\rm de}$).

While numerous statistics of cosmic structures have been used to understand
the evolution and structure formation of our universe, the exact definition of
cosmic structures remains unclear. This is mainly because the matter
distribution on large scale is continuous, and therefore, there exists no
specific discrete boundary for each structure. Also, the membership for a
certain structure might change if one considers other properties than just
position, such as dynamics, mass, and so on (e.g., Serra & Diaferio, 2013;
Gifford et al., 2013). Due to this ambiguity, numerous clustering algorithms
have been proposed and used in the astronomical community. For example, Knebe
et al. (2011) compared the various properties of dark matter (DM) halos and
subhalos found by 17 different halo-finding algorithms run on the same
cosmological $N$-body simulation and found them to be in close agreement.

Galaxy clustering algorithms have played a crucial role in identifying galaxy
clusters or super-clusters and exploring the large-scale structure of the
universe, including filamentary structures. The most commonly utilized
algorithms in astronomical research are the friends-of-friends (FoF; Davis et
al. (1985)) and the minimum spanning tree (MST; Borůvka (1926)) algorithms,
which were introduced in the 1980s and have since become the standard for
galaxy clustering. However, more recent algorithms, such as the halo-based
algorithm (Yang et al., 2007) and redMapper (Rykoff et al., 2014), have also
been developed to address some of the limitations of the FoF and MST
algorithms and are more directly applicable to observational data.

In recent years, with the rapid development of machine learning (ML)
technology, clustering algorithms such as DBSCAN (Density-based Spatial
Clustering of Applications with Noise; Ester et al., 1996) have also been
applied to galaxy clustering. These clustering softwares, including the ML
based one, show comparable performance in galaxy clustering and produce
consistent clustering results. However, the results of clustering do not
always represent the clusters that the human eye can find. There are cases
where a distribution clearly contains a cluster but it is not recognized as
such by clustering algorithms, and there are cases where it is clearly divided
into two clusters, visually, but appears as a single lump to the software.

As such, we have explored the possibility of developing an algorithm that
creates galaxy clustering (or, networks) in a way that more closely resembles
how the human eye and brain identify patterns. One approach we considered was
to adapt jet-finding software used in high-energy particle physics research,
with a particular focus on the MulGuisin (MGS) algorithm as a potentially
suitable software for galaxy clustering. MulGuisin (물귀신) is a Korean word for
a ghost that lives in water and is a figure that often appears in old Korean
stories. The MGS algorithm started with the idea that the ghosts hiding in the
water could be found in the order of height by simply draining the water from
the lake.

Initially, we copied the MGS software from the ATLAS (_A Toroidal LHC
Apparatus_) Jet-Finding library released in the early 1990s and developed it
into a 3D galaxy network finding algorithm. We then made several sample galaxy
distributions to check the performance of MGS and compared the results to
those produced by other standard clustering algorithms such as FoF, MST, and
DBSCAN. As a result, it was found that MGS has several important
characteristics that other algorithms cannot demonstrate. The most important
of these features is that as the linking length parameter for clustering
increases, percolation occurs in all other algorithms, whereas in MGS, such a
phenomenon does not occur even when a very large linking length is used.
Therefore, MGS can be useful for finding large overdensity regions that are
much broader than typical galaxy clusters by using very large linking lengths,
which are not possible with other algorithms.

The networks identified by MGS are not necessarily all distinct,
gravitationally bound collections of galaxies. Rather MGS is a general method
for identifying network structures in a spatial point set. In this work we
look at the performance of MGS at identifying clustered galaxies, compared
with other algorithm, e.g. FoF, etc. In subsequent work we will look in more
detail at the physical properties of the MGS networks, e.g. halo mass
function, and their usefulness in cosmological analysis.

The structure of this paper is as follows. In Section 2, we introduce the
MulGuisin clustering algorithm, as well as FoF, MST, and DBSCAN as benchmark
clustering algorithms for comparison. In Section 3, we describe both
controlled random data and realistic galaxy distribution data that we will use
for the performance test. We apply the above four clustering algorithms to the
data and compare their performances in Section 4, and we summarize our results
in Section 5.

##  2 Methods

In this section, we first introduce our MulGuisin clustering algorithm and
introduce two clustering tools that are widely used in the field of astronomy
and a newly developed clustering program through machine learning.

###  2.1 MulGuisin network finding algorithm

Figure 1: Schematic flow chart to describe how MulGuisin (MGS) algorithm works

The MulGuisin (MGS) network finding algorithm was first introduced as a jet
finder for the _Large Hadron Collider_ (LHC) physics in the ATLAS
Collaboration (Bosman et al., 1998). The algorithm is neither a variant of the
conventional cone jet finding algorithm nor a variant of the $k_{T}$
clustering algorithm that is used in various collider experiments as the
standard tools for finding jets. Although it has shown some improvements in
jet reconstruction performance, such as optimized jet direction and jet energy
resolution, but has not been used as a standard jet-finding tool for LHC
experiments.

Fig. 1 shows how the MGS algorithm works. The MGS algorithm first finds the
point with the highest local number density from the input data and designates
it as a network seed. Then it finds the point with the second highest local
number density and decides whether the point should belong to the first
network or become a seed of a stand a new network. This decision is made by
checking how close the test point is to any neighboring networks, for which we
introduce a parameter called the linking length ($\ell_{\rm Link}$). That is,
if the distance between the test point and the closest point in the network is
less than this parameter, the test point is attached to the network,
otherwise, it becomes the seed of a new network. This parameter can be called
resolution, which governs the separation of networks. The algorithm then finds
the next massive point and repeats the above process until there are no more
points left to test. At this stage, all points are assigned into networks or
remain as isolated points.

Figure 2: Diagram showing how the MulGuisin (MGS) algorithm works to identify
networks. Each gray circle represents a galaxy, and the size of the circle
denotes its local number density, with the number specifying the galaxies
ranking in descending order of density. In this specific example, 21 galaxies
are grouped into 3 networks and 2 isolated galaxies.

Fig. 2 is an illustration to explain how the MGS algorithm creates galaxy
networks. In the figure, the points are sorted in order according to their
local number density which is estimated using a Voronoi tessellation. And
according to this order, they become new network seeds or become attached to
existing networks. After going through the process, a network forms a tree-
like structure that is sequentially connected according to the order of local
number density. The points in a network then form branches and nodes, and from
the characteristic structure of such tree shape, one may able to study the
topology of the galaxy network.

###  2.2 Benchmark Algorithms

In order to compare the performance of the MGS as a new clustering tool with
the typical clustering algorithms, we select three benchmark algorithms,
mostly based on their popularity in the astronomical community, mathematical
clarity, and versatility. They are the Friends-of-Friends (FoF), Minimum
Spanning Tree (MST), and the Density-Based Spatial Clustering of Applications
with Noise (DBSCAN). Here we briefly introduce each package and describe how
they make networks.111Note that running our MGS algorithm from scratch may
take a longer time than the above benchmark algorithms, especially when the
number of data points is large. We found that most of the MGS calculation
time, for a large number of data points, is taken in constructing the Voronoi
tesselation and calculating the local density for each point. If we separate
MGS into a density calculation and a tree building part, we found that the
tree building takes a similar time to the benchmark algorithms.

####  2.2.1 Friends-of-Friends (FoF)

The friends-of-friends (FoF) algorithm is a commonly used technique for
identifying networks in astrophysical data (Huchra & Geller, 1982; Tago et
al., 2008; Duarte & Mamon, 2014; Tempel et al., 2016). This algorithm has a
single free parameter, the linking length ($\ell_{\rm FoF}$), which determines
the distance threshold for linking two data points. Points that are within
this distance of each other are considered to be connected, and all connected
points are grouped together into a single cluster.

One limitation of the FoF algorithm is that it can be difficult to choose an
appropriate linking length. Different values of this parameter can result in
networks of different shapes or numbers, making it challenging to determine
the optimal value (Tago et al., 2008).222Note that the appropriate choice
linking length for identifying DM halos from the DM particles in the $N$-body
simulations is well known ($\ell_{\rm FoF}\simeq 0.2\langle d_{\rm
particle}\rangle$) (More et al., 2011). However, the optimal choice of linking
length in general clustering problems is not well known. In this study, we use
the Halotools implementation of the FoF algorithm (Hearin et al., 2017) to
identify networks in our datasets by applying various $\ell_{\rm FoF}$. Unless
otherwise noted, we assume all FoF groups containing two or more members as
networks.

####  2.2.2 Minimum Spanning Tree (MST)

Galaxy data can be represented as a graph, with each galaxy represented as a
node and the distance between two galaxies represented as an edge. The minimum
spanning tree (MST) algorithm is a method for constructing a unique network
from this data by connecting all nodes with minimum edges. Unlike other
clustering algorithms, the MST does not require the use of a free parameter
such as a linking length to construct the entire network. However, the MST
connects all nodes and may not produce networks with shapes that accurately
reflect those of the original networks.

Nevertheless, MST has been used in cosmology to study the large-scale
structure of the universe (Barrow et al., 1985; Krzewina & Saslaw, 1996;
Naidoo et al., 2020). In this study, we use the MiSTree package (Naidoo, 2019)
to construct MSTs from our galaxy data. Then, we find networks from the single
MST tree by cutting nodes longer than the linking length ($\ell_{\rm MST}$).
Similar to the FoF case, we apply various values of $\ell_{\rm MST}$ and
assume all tree segments containing two or more members as networks.

####  2.2.3 Density-based Spatial Clustering of Applications with Noise
(DBSCAN)

The use of machine learning (ML) techniques is widespread in astronomy, as
they enable the identification of patterns in data using algorithms. ML
algorithms can be classified based on the type of data they are applied to,
and one type, called unsupervised ML, is used with unlabeled data. Clustering
algorithms, a subcategory of unsupervised ML algorithms, are used to group
together data points with similar properties. One popular clustering algorithm
is DBSCAN (density-based spatial clustering of applications with noise), which
has been applied in a variety of contexts (Ester et al., 1996; Sander et al.,
2017).

DBSCAN is a density-based clustering algorithm that groups together data
points based on their local density. In this algorithm, each network is
identified by defining its core, which consists of high-density points within
a certain distance. The definition of core requires two free parameters,
min_samples and eps, which determine the minimum number of neighbors a point
must have within a given radius in order to be considered as the core. Then,
other points that are directly reachable from some core points within eps are
also considered part of the network, while other points are labeled as noise.

In this study, we use the scikit-learn package (Pedregosa et al., 2011) to
implement the DBSCAN algorithm and identify networks in our data by applying
various eps (or, the “linking length” in DBSCAN ($\ell_{\rm DBSCAN}$)). Unless
otherwise noted, we assume min_samples $=3$.

###  2.3 A Simple 2D Toy Model Test

To see how the shape of the networks generated by the MGS algorithm differs
from the results of other clustering algorithms, we created simple simulation
data and compared the results. We first assume that there are 5 networks in 2D
space, and consider the case where each network contains 50 galaxies equally.
The galaxies are distributed within the network according to a Gaussian
function with a width of $10\mathinner{h^{-1}{\rm Mpc}}$. The coordinates of
the two-dimensional space span from 0 to $100\mathinner{h^{-1}{\rm Mpc}}$ on
both the X- and Y-axes, and the position of each network is set to have three
different distributions, from far away from each other to all close together,
as shown in Fig. 3 column (a) in rows (1), (2) and (3).

Figure 3: A simple 2D toy model test of the MGS algorithm by comparing it with
the MST algorithm. (a) Input distributions of 5 networks with different
degrees of separation from each other ((1)–(3)). Background color denotes the
galaxy number distribution we used for generating the galaxies. (b) Networks
found by MGS and their tree structures. (c) Networks found by MST and their
tree structures. We used a linking length (LL) of $10\mathinner{h^{-1}{\rm
Mpc}}$ for both algorithms and we draw the size of linking length in the top
middle figure.

As shown in Fig. 3, both the MGS and MST algorithms correctly find 5 networks
when the distances among the networks are sufficiently far apart. However,
when the networks get closer together, MST cannot differentiate between the
networks and starts recognizing them as one big network. Even for the cases
where networks are attached to each other as shown in Fig. 3 (3), MGS still
recognizes four among five true networks like the human eyes can distinguish
each network, whereas MST recognizes 4 adjacent networks as one huge network.
These differences can create serious differences in results when studying the
number and mass distributions of networks.

In the next section, we will compare the performance of MGS and other
algorithms with more realistic 3D data.

##  3 Data

Data set |  Description
---|---
D1-LD |  50 randomly positioned networks, each of which contains 100 galaxies randomly spread by the 3D Gaussian distribution with standard deviation $\sigma=1\mathinner{h^{-1}{\rm Mpc}}$. The total number of galaxies is 5,000.
D1-HD |  Same as D1-LD, but with the greater standard deviation $\sigma=10\mathinner{h^{-1}{\rm Mpc}}$.
D2-NA |  Same as D1-HD, but the number of galaxies in each network follows an exponential random distribution. The total number of galaxies is 7,041.
D2-LA |  Same as D2-NA, but adding uniformly randomly distributed noisy galaxies to the entire box to increase the total galaxy number density 1.5 times of D2-NA.The total number of galaxies is 12,041.
D2-HA |  Same as D2-LA, but adding more noisy galaxies so that the total galaxy number is twice D2-NA. The total number of galaxies is 17,041.
D3-HOD |  500 randomly positioned networks with a mass distribution similar to the Press-Schechter mass function. The number of galaxies for each network follows HODa for massive halos ($M\geq 10^{13}\mathinner{h^{-1}{\rm M}_{\odot}}$). The galaxies are spread by the NFW profile with the concentration parameter to 10. The total number of galaxies is 65,460.
D4-SDSS |  Volume-limited sample of the KIAS-VAGCb with absolute $r$-band magnitude $\mathcal{M}_{r}-5\log h<-20$.
D4-HR4 |  Four lightcone data of mock galaxy catalogs from the Horizon Run 4 simulationc with a similar condition to D4-SDSS.

a Kravtsov et al. (2004). b Choi et al. (2010a). c Kim et al. (2015); Hong et
al. (2016).

Table 1: Name and description of galaxy data sets that we use in this
analysis. The box size of all controlled data (D1–D3) is
$(200\mathinner{h^{-1}{\rm Mpc}})^{3}$.

Our final goal is to use the MGS algorithm for the galaxy clustering or other
large-scale structure finding for the cosmology study. However, since some
inconsistencies exist between various clustering algorithms for finding
networks or other large-scale structures (e.g., see Knebe et al., 2011, and
references therein), we cannot compare the MGS networks found in the realistic
data with their “truth”.

Therefore, we apply two types of data sets in this section to compare the
performance between MGS and other benchmark algorithms. The first sets, called
the “controlled random data” (D1–D3), are those that we design all properties
of networks, including their positions and member galaxy distributions. Since
we already know the true information of each network, we can test which
algorithms predict the true networks better in which conditions. The next
sets, called the “realistic data” (D4), are the observational and simulation
data sets of galaxies around $z\simeq 0$, and we focus on comparing the
properties of predicted networks in each algorithm. Table 1 summarizes the
data sets we use in this work.

###  3.1 Controlled Random Data (D1–D3)

Figure 4: Three-dimensional galaxy distributions of controlled data sets used
in this paper. Top: D1-LD(left), D1-HD(middle), and D3-HOD(right). Bottom:
D2-NA(left), D2-LA(middle), and D2-HA(right), with noisy additional galaxies
shown as yellow dots. See Table 1 for details.

####  3.1.1 Different Spatial Dispersion (D1)

We use controlled, simulated data to evaluate the performance of the MGS
algorithm in comparison to other clustering algorithms. These data are
generated randomly and allow us to control the shape and distribution of
networks to test the algorithms under different conditions. The first set of
data consists of 100 galaxies per network and 50 networks within a
3-dimensional cubic volume of space with a side length
$200\mathinner{h^{-1}{\rm Mpc}}$. The network center positions are chosen
randomly, and the galaxies in each network are distributed according to a
Gaussian distribution with a variable standard deviation ($\sigma$) that
controls the spatial dispersion. The D1-LD data has a low spatial dispersion
($\sigma=1\mathinner{h^{-1}{\rm Mpc}}$), leading to well-separated networks,
while the D1-HD data has a higher spatial dispersion
($\sigma=10\mathinner{h^{-1}{\rm Mpc}}$), resulting in networks that are
closer together. Upper left and middle panels of Fig. 4 show the distribution
of galaxies in the D1-LD and D1-HD data sets.

####  3.1.2 Additional Noisy Galaxies (D2)

We generate additional controlled data sets that are similar to D1 but with
slightly different characteristics. We again place 50 network centers at the
same positions as D1, but this time we use an exponential distribution to
generate a variable number of galaxies for each network:

$P(N_{\rm gal})=\left\\{\begin{array}[]{ll}\displaystyle{\frac{1}{\Delta
N}\exp\left[-\frac{N_{\rm gal}-N_{0}}{\Delta N}\right]}&\textrm{ if }N_{\rm
gal}>N_{0}\\\ 0&\textrm{ otherwise}\\\ \end{array}\right.~{}.$ (1)

Here, we set $N_{0}$ and $\Delta N$ as 50 and 100, respectively, so that the
minimum number of galaxies per network and the total number of galaxies
roughly match with D1. The galaxies are spatially distributed according to a
Gaussian function centred on the network’s centre with a standard deviation of
$\sigma=10\mathinner{h^{-1}{\rm Mpc}}$, as was done for D1-HD, and we call
this new controlled data D2-NA. The lower left panel of Fig. 4 shows the
distribution of galaxies in the D2-NA data. The total number of galaxies in
this data set is 7,041.

In addition to D2-NA, we introduce two more data sets that are created by
adding unclustered galaxies, which are sampled uniformly in the entire box. We
add these ‘noisy’ galaxies so as to test how the algorithms are affected by
the background density. The lower middle panel of Fig. 4 shows the D2-LA data,
where we add 5,000 galaxies (yellow dots) to increase the galaxy number
density by 1.5 times compared to D2-NA. On the other hand, the lower right
panel shows the D2-HA data, where we add 10,000 galaxies to make the total
galaxy number density twice that of D2-NA.

####  3.1.3 HOD-based Mock Galaxies (D3-HOD)

We generate a third set of controlled, simulated data to create a more complex
environment for testing the performance of the MGS algorithm. We use an
analytic formula to model the distribution of galaxies in this data set.
First, we create 500 network center positions by sampling uniform random
distribution within a $(200\mathinner{h^{-1}{\rm Mpc}})^{3}$ box. Then, we
obtain the normalized version of Press-Schechter halo mass function at $z=0$
(Press & Schechter, 1974), with a concordance $\Lambda$CDM cosmology to the
Planck 2015 data (Planck Collaboration et al., 2016), for massive halos
$M_{\rm halo}>10^{13}\mathinner{h^{-1}{\rm M}_{\odot}}$ using the Colossus
package (Diemer, 2018). We then obtain masses for each of the 500 networks by
randomly sampling for the mass function.333Note that neither the positions nor
the mass distribution of networks in D3-HOD follows the estimation from the
standard cosmology. However, here we focus only on providing complex
environments, and therefore, such differences do not affect our motivation.
See Section 3.2 for realistic data sets instead.

We then use this information to generate a distribution of the number of
galaxies using a halo occupation distribution (HOD) model. The mean halo
occupation is typically assumed to follow a power law at massivehalo masses
(Berlind & Weinberg, 2002; Kravtsov et al., 2004):

$N_{\mathrm{avg}}(M_{\rm
network})=\left\\{\begin{array}[]{ll}\displaystyle{\left(\frac{M_{\rm
network}}{M_{1}}\right)^{\alpha}}&\textrm{if }M_{\rm halo}>M_{\mathrm{min}}\\\
0&\textrm{otherwise}\\\ \end{array}\right.~{},$ (2)

where $\alpha$, $M_{\rm min}$, and $M_{1}$ correspond to the power-law index,
cutoff halo mass where halo cannot contain galaxies, and the mass scale
containing a single galaxy at the given condition of galaxy sample. Here, we
use $\alpha=0.87$ and $M_{\mathrm{min}}=10^{13}\mathinner{h^{-1}{\rm
M}_{\odot}}$ by following Kravtsov et al. (2004). We set
$M_{1}=10^{11}\mathinner{h^{-1}{\rm M}_{\odot}}$ so that the minimum number of
galaxies for each network is set as 50. Also, for simplicity, we calculate the
actual number of galaxies at each network by applying the ceiling to $N_{\rm
avg}$.

Next, we use the Colossus package to create an Navarro-Frenk-White (NFW)
profile (Navarro et al., 1996)

$\rho(x)=\frac{M_{\rm network}}{4\pi R_{\rm
vir}^{3}}\left[\left\\{\ln(1+c_{\rm s})-\frac{c_{\rm s}}{1+c_{\rm
s}}\right\\}x(x+c_{\rm s}^{-1})^{2}\right]^{-1}~{},$ (3)

where $x\equiv r/R_{\rm vir}$. The concentration parameter for the NFW profile
$c_{\rm s}$ is fixed as 10, and the virial radii $R_{\rm vir}$ is determined
by the network mass accordingly. We then randomly distribute the galaxies
according to this profile, and the resulting data set consists of 65,460
galaxies. The upper right panel of Fig. 4 shows the distribution of galaxies
in D3-HOD.

###  3.2 Realistic Data: SDSS & Horizon Run 4 (D4)

In the previous subsection, we described a set of controlled data catalogues
for which we can carefully control the properties of the networks. Such data
are useful for testing the performance of MGS over other benchmark algorithms
by comparing the properties of identified networks with the input truth.
However, the true distribution of galaxies in the universe differs from these
controlled random data in the following ways. First, unlike those in the
controlled random data with low noise levels, the boundaries of clusters in
the universe are often not clearly defined (e.g., Serra & Diaferio, 2013;
Gifford et al., 2013). Also, the spatial distribution of galaxies in each
cluster may not follow spherical symmetry (e.g., Limousin et al., 2013, for a
good review). Furthermore, the redshift-space distortion elongates spherical
clusters in real space, which may require that we separate linking lengths
between the radial and tangential directions (Farrens et al., 2011; Tempel et
al., 2016).

Therefore, it is necessary to adopt a realistic galaxy distribution for a fair
performance test of the MGS algorithm. However, unlike for the case of the
controlled random data where we know the answer, we may only study the
difference between the network properties from the MGS and other benchmark
algorithms. Here we use observational data and four corresponding sets of mock
simulation data — the volume-limited KIAS-Value Added Galaxy Catalog (KIAS-
VAGC) of the _Sloan Digital Sky Survey_ (SDSS) Main Galaxy Sample with
$r$-band absolute magnitude $\mathcal{M}_{r}-5\log h<-20$ (Choi et al., 2010a)
and the lightcone mock galaxy samples from the Horizon Run 4 simulation (Kim
et al., 2015; Hong et al., 2016).

####  3.2.1 Volume-limited KIAS-VAGC (D4-SDSS)

The KIAS Value-Added Galaxy Catalog (KIAS-VAGC; Choi et al., 2010a) is an
upgraded version of the New York University Value-Added Galaxy Catalog (NYU-
VAGC; Blanton et al., 2005), which is part of the _Sloan Digital Sky Survey_
(SDSS) Data Release 7 (Abazajian et al., 2009), by adding some missing
redshifts to improve spectroscopic completeness.444The completeness of
spectroscopic redshift for our sample (D4-SDSS) is 94.2%. The remaining 5.8%
had suffered the fiber collision(see the main text), and the spectroscopic
redshifts of their nearest neighbor galaxies are assigned to them instead.
This catalog has been widely used in numerous studies, including cosmic voids
statistics (Pan et al., 2012; Hoyle et al., 2012), largest structures of
universe (Park et al., 2012), fraction of barred galaxies (Lee et al., 2012a),
and the properties of active galactic nuclei (AGN; Lee et al., 2012b; Hwang et
al., 2012; Bae & Woo, 2014).

Figure 5: Selection of the volume-limited sample of the KIAS Value-Added
Galaxy Catalog (KIAS-VAGS) used in this study (red boxes). Left: Volume-
limited selection in the redshift vs. absolute $r$-band magnitude plane with
$\mathcal{M}_{r}-5\log h<-20$. Middle & Right: Sky selection shown in both
SDSS Survey coordinates $(\eta,\lambda)$ (middle) and (RA, DEC) coordinates
(right).

Most of the KIAS-VAGC galaxies were observed with the apparent $r$-band
magnitude limit $r=17.6$. It means that, in terms of absolute magnitude, the
catalog contains less bright galaxies at lower redshifts, while only very
bright galaxies could be seen at higher redshifts. Therefore, for a fair
comparison between galaxies over a wide redshift range, we apply a “volume-
limited” selection by selecting galaxies brighter than a certain absolute
$r$-band magnitude (Choi et al., 2010b). Here, we use $\mathcal{M}_{r}-5\log
h<-20$. By combining with the given apparent $r$-band magnitude limit, such
absolute magnitude cutoff naturally provides the upper redshift bound of our
volume-limited sample ($z<0.107$; left panel of Fig. 5). We also apply the
lower redshift bound $z>0.02$, by considering the incompleteness of the galaxy
sample below the given redshift.

In addition to the volume-limited selection in the redshift-magnitude plane,
we also apply a sky selection for simplification. Specifically, instead of
using spherical coordinates with right ascension (RA) and declination (DEC)
$(\alpha,\delta)$, we select galaxies with the SDSS Survey coordinate ($\eta$,
$\lambda$), defined by

$\displaystyle\eta$ $\displaystyle\equiv$
$\displaystyle\tan^{-1}\left[\frac{\sin\delta}{\cos\delta\cos(\alpha-\alpha_{0})}\right]-\delta_{0}$
(4) $\displaystyle\lambda$ $\displaystyle\equiv$
$\displaystyle\sin^{-1}\left[\cos\delta\sin(\alpha-\alpha_{0})\right]\,.$ (5)

Here, $(\alpha_{0},\delta_{0})=(185^{\circ},32.5^{\circ})$ is the RA/DEC sky
position where $(\eta,\lambda)=(0,0)$. We select the rectangular sky region
$-33.5^{\circ}<\eta<36.5^{\circ}$ and $-48^{\circ}<\lambda<51^{\circ}$, in
order to maximize the sky area with a simple geometry, and to avoid issues
arising from a complicated boundary (middle and right panels of Fig. 5).

####  3.2.2 Horizon Run 4 (D4-HR4)

The Horizon Run 4 simulation (HR4; Kim et al., 2015) is an extremely large
cosmological $N$-body simulation that uses $6,300^{3}$ DM particles within a
periodic cube with a comoving volume $V=(3.15~{}h^{-1}{\rm cGpc})^{3}$. It
assumes a vanilla $\Lambda$CDM cosmological model in concordance with the
_Wilkinson Microwave Anisotropy Probe_ (WMAP) 5th-year result (Dunkley et al.,
2009). Among 2,001 timesteps between $z=100$ to 0, 75 coarse timesteps with
mean time difference $\Delta t=0.18~{}{\rm Gyr}$ are chosen between $z=12$ to
0 to build a merging tree of FoF halos. The FoF linking length is $0.2$ times
the particle mean separation, and we identify halos only whose mass is greater
than $M_{\rm halo}^{\min}=2.7\times 10^{11}\mathinner{h^{-1}{\rm M}_{\odot}}$.

The mock galaxies are then produced by so-called the most bound halo particle
(MBP)-galaxy abundance matching method (Hong et al., 2016). We find MBPs for
all halos in the merging tree and adopt their positions and peculiar
velocities as those of corresponding mock galaxies. The “mass” of mock
galaxies, which is used as a proxy of stellar mass or luminosity, is defined
as the mass of their hosting halos. For satellite halos, we identify their
MBPs at the timestep just before the infall event and trace them until they
are totally absorbed toward their central halo by tidal disruption. For
estimating the tidal disruption timescale $t_{\rm merge}$), we adopt a
modified model of Jiang et al. (2008),

$\frac{t_{\rm merge}}{t_{\rm
dyn}}=\frac{(0.94\epsilon^{0.60}+0.60)/0.86}{\ln[1+(M_{\rm host}/M_{\rm
sat})]}\left(\frac{M_{\rm host}}{M_{\rm sat}}\right)^{\alpha}~{},$ (6)

where $\epsilon$, $M_{\rm host}$, $M_{\rm sat}$, $t_{\rm dyn}$ are the
circularity of the satellite’s orbit, the mass of central and satellite halos,
and the orbital period of virialized objects, respectively. We adopt
$\alpha=1.5$ for a better match of the galaxy two-point correlation function
(2pCF) at scales less than $1\mathinner{h^{-1}{\rm Mpc}}$ at a given spatial
resolution of the HR4 (Zehavi et al., 2011; Park et al., 2019). Then the mass
of survived satellite galaxies is defined as the mass of their hosting halos
just before the infall.

After producing snapshot mock galaxy catalogs for coarse timesteps, we then
produce lightcone mock galaxy catalogs up to $z=1.5$. The all-sky lightcone DM
particle data of the HR4 were created during the simulation by stacking the
comoving shells at the corresponding redshifts. Then, we compare the IDs of
the galaxy MBPs at each coarse timestep snapshots and those of DM particles at
the lightcone data with the coarse comoving shells. If the MBP ID of a given
mock galaxy matches that of a particle in the lightcone data, we assign a
galaxy in the lightcone data. Here, we adopt the position and peculiar
velocity from the particle at the lightcone data, while the galaxy “mass”
comes from the mock galaxy at the nearest snapshot.

After creating the all-sky lightcone mock galaxy catalog, we cut it in a
similar way to the volume-limited KIAS-VAGC sample. First, we apply the
redshift space distortion (RSD) for each mock galaxy for a fair comparison
with observation, by using real-space positions and peculiar velocities. Then,
we apply the same redshift range $0.02<z<0.107$ and set the lower bound of
galaxy “mass,” so that the galaxy number density of the HR4 lightcone data is
identical to that of KIAS-VAGC. After that, we create four non-overlapping
subsets from it with the same angular geometry as our SDSS Survey coordinate
selection.

During the analysis, we found that the fiber collision in the fiber-fed
spectroscopic observations affects various clustering statistics (Zehavi et
al., 2002; Guo et al., 2012; Reid et al., 2014; Tonegawa et al., 2020).
Therefore, for a fair comparison, our HR4 mock galaxy catalogs also need to
follow the same fiber collision condition as the KIAS-VAGC. To do so, we
select pairs of mock galaxies whose angular distance is less than 55
arcseconds and keep only one from each pair by random selection. Because SDSS
observations were partially overlapping, some close-pairs have both redshifts.
In order to reflect this, we only fiber-collide 60% of the close pairs.

##  4 Results

We evaluate the performance of the MGS algorithm and other relevant methods
using three controlled datasets and observational data. Four distinct
algorithms are executed with varying linking-lengths to determine the number
of networks. However, assessing these algorithms solely based on the number of
networks identified is insufficient. Even if an algorithm successfully
recovers the correct network count, it remains uncertain whether the
identified networks genuinely correspond to the original networks.

To account for this property and ensure that the estimated networks from the
algorithms align with the original networks, we employ the Adjusted Rand Index
(ARI; Morey & Agresti, 1984; Hubert & Arabie, 1985), a modification of the
Rand Index (RI; Rand, 1971). First, RI measures the similarity between two
clustering outcomes by counting pairs of elements under specific conditions
and is defined as

$\mathrm{RI}=\frac{\mathrm{TP}+\mathrm{TN}}{\mathrm{TP}+\mathrm{FP}+\mathrm{FN}+\mathrm{FP}}.$
(7)

Here, TP (true positives) represents the number of pairs where the algorithm
accurately identifies original network members. TN (true negatives) denotes
the number of pairs where the algorithm correctly distinguishes non-members.
FP (false positives) stands for the number of pairs where the algorithm
incorrectly classifies non-members as members. FN (false negatives) signifies
the number of pairs where the algorithm wrongly classifies members as non-
members.

The ARI is then defined as

$\mathrm{ARI}=\frac{\mathrm{RI}-\mathbb{E}[\mathrm{RI}]}{1-\mathbb{E}[\mathrm{RI}]}.$
(8)

Here, $\mathbb{E}[\mathrm{RI}]$ is an expectation value of the RI upon
permutation of clustering indices. An ARI value close or less than zero
indicates that the algorithm is unable to recover the original network, while
a value close or equal to one signifies efficient recovery of the original
network. We assess the performance of the algorithms using ARI as a function
of linking length. Here, we set the minimum value of ARI as zero.

We use the three controlled data sets and observation data to evaluate the
performance of the MGS algorithm and compare it to other clustering
algorithms. We run each of the four algorithms with different values of the
linking length and calculate the ARI of results identified by each algorithm.
We repeat this process for a range of linking lengths and analyze the trends
in the number of networks identified by each algorithm. This allows us to
assess the sensitivity of the algorithms to the choice of linking length and
to compare their performance in identifying networks in the different data
sets.

###  4.1 Results with Controlled Data

Figure 6: The Adjusted Rand Index (ARI) as a function of linking length for
D1-LD (top panel) and D1-HD (bottom). Each of the 4 clustering algorithms are
indicated using different colors and symbols. Note that MST and DBSCAN show
considerable overlap.

Fig. 6 shows the ARI obtained by each algorithm as a function of the linking
length for the controlled data set 1. The top panel shows the results for the
D1-LD data, which consists of well-separated networks. All four algorithms
perform well in identifying the networks, but the MGS algorithm stands out for
its ability to accurately identify the correct number of networks. In
particular, for large linking lengths, the FoF and DBSCAN algorithms identify
fewer than 50 networks, resulting in a lower ARI, because they connect
neighboring networks and merge them into a single network.

The bottom panel of Fig. 6 shows the results for the D1-HD data, which has a
higher level of spatial dispersion and some networks that are close to each
other. For smaller linking lengths, the algorithms identify many more networks
than the true input value of 50. If the linkage length is too small the
algorithm cannot connect all the galaxies in each network, and only a few
galaxies become a new network. That is, one network is artificially divided
into many individual segments, and consequently the algorithms have very low
ARI value. As the linking length increases, the behavior of the algorithms
becomes more distinct. The MGS algorithm continues to accurately identify the
correct number of networks, while the other algorithms identify fewer networks
due to the merging of originally separate networks. This is the opposite
reason for why the ARI is low at larger linking lengths as compared lower
linkage lengths. The MGS algorithm is able to track the structure of the
networks and identify their boundaries, leading to more accurate results in
this type of data.

Figure 7: Same as Fig. 6, but with D2-NA(top), D2-LA(middle), and
D2-HA(bottom). The vertical dashed line is mean-separation of data ($\langle
d\rangle$). Since each data set has a different overall number density. Figure
8: 3D distribution of clustering results from MGS and other benchmark
algorithms in D2-LA with linking length $\ell=11\mathinner{h^{-1}{\rm Mpc}}$,
which is longer than the mean-separation $\langle
d\rangle=8.73\mathinner{h^{-1}{\rm Mpc}}$. Color indicates the network
membership. MGS finds 49 networks among 50 true networks, while other
algorithms connect most of galaxies and finally make a giant network (blue
color).

Fig. 7 shows the results for controlled data set 2, which includes the D2-NA
data with no additional galaxies and the D2-LA and D2-HA data with additional
galaxies. The top panel shows the ARI identified by each algorithm for the
D2-NA data. When this data was generated, the minimum number of galaxies per
network was set to 50. In the region of small linking lengths, all algorithms
identify fewer than 50 networks and obtained low ARI because the linking
length is too small to connect the galaxies in the networks. As a result, the
networks identified by the algorithms have fewer than 50 member galaxies, and
are therefore not considered as true networks. For larger linking lengths,
particularly those larger than 5, the difference between the MGS algorithm and
the other algorithms becomes more pronounced. The MGS algorithm continues to
accurately identify the correct number of networks, while the other algorithms
identify fewer networks due to the merging of originally separate networks.

The behavior of the algorithms with additional galaxies is even more distinct.
The middle panel of Fig. 7 shows the results for the D2-LA data, where the
other three algorithms identify only a single network for very large linking
lengths. As the linking length increases, the algorithms merge several
networks into a single giant network, resulting in a significantly lower
number of networks than the original data and consequently a low ARI value.
This rapid increment of a single giant network is called “percolation,” and it
is known to occur at linking length similar to the mean-separation
($\ell\simeq\langle d\rangle$) for the ideal random Poisson graph (Dall &
Christensen, 2002). Fig. 7 clearly shows that such percolation occurs at
$\ell\simeq\langle d\rangle$ for all three benchmark algorithms.

Note that, however, FoF has the shortest linking length at percolation, while
both MST and DBSCAN share a similar value of linking length at percolation.
This is because FoF does not have an additional consideration for limiting the
network boundary that exists in the other two algorithms (minimize the number
of edges in MST, and core definition in DBSCAN). Fig. 8 shows the 3D
distributions of clustering results from various algorithms at linking length
$\ell=11\mathinner{h^{-1}{\rm Mpc}}$, which is longer than the mean separation
$\langle d\rangle=8.73\mathinner{h^{-1}{\rm Mpc}}$. As expected, three
benchmark algorithms show percolation (blue color), while our MGS algorithm
successfully reconstructs most of the true networks. Note that only one giant
network is found in the FoF algorithm, while both MST and DBSCAN have two
additional small networks (green and orange colors).

The behavior of the MGS algorithm for the D2-HA data is slightly different. In
the region of small linking lengths, the MGS algorithm accurately identifies
the correct number of networks. However, for larger linking lengths,
particularly $\ell>13\mathinner{h^{-1}{\rm Mpc}}$, the MGS algorithm
identifies additional networks that were not present in the original data
(corresponding to a low ARI value). These “fake” networks are not true
networks and are not representative of the underlying structure of the data.
This behavior highlights the ability of the MGS algorithm to identify networks
in data with a complex distribution of galaxies but also underscores the
importance of choosing an appropriate linking length to avoid identifying
false networks.

Figure 9: Same as Figs. 6–7, but with D3-HOD.

Fig. 9 shows the number of networks for controlled 3 data with a more complex
environment than D1–D2. At $\ell\ga\langle d\rangle/2\approx
3\mathinner{h^{-1}{\rm Mpc}}$, the number of networks and the ARI value using
FoF and DBSCAN decreases as the linking length increases, resulting in the
percolation at $\ell\ga\langle d\rangle$. The MST shows a flat curve when the
linking length is larger than $\sim 8\mathinner{h^{-1}{\rm Mpc}}$. This is
because MST connects all galaxies with minimal edge first, and then we cut off
the links with linking length. Therefore, if there were no links longer than
$8\mathinner{h^{-1}{\rm Mpc}}$ in the original tree, then cutting the links
with any longer linking length than $8\mathinner{h^{-1}{\rm Mpc}}$ would not
change the result. So, the ARI value using MST shows a constant value.

In contrast, the ARI value and also the number of networks identified by the
MGS algorithm slowly decreases as the linking length increases. This is
because the networks in this data set are close to each other and are easily
merged by the algorithm for large linking lengths. However, the MGS algorithm
is able to identify networks based on density, which allows it to retain the
structure of the networks even for large linking lengths. This is the main
advantage of the MGS algorithm compared to the other three algorithms, which
are not able to accurately identify networks in complex data sets.

###  4.2 Results with Observational and Cosmological Simulation Data

Figure 10: The number of networks as a function of linking length for D4-SDSS
(thick lines) and D4-HR4. For D4-HR4, the average values and the ranges
between minimums and the maximums of 4 data samples are drawn as thin lines
and error bars. Each panel shows the results from the different choices of the
minimum number of member galaxies to identify networks.

Fig. 10 shows the results of the four algorithms applied to both KIAS-VAGC
observational data and four sets of HR4 lightcone data. In the previous
subsections, we calculated the ARI in order to evaluate the algorithms. One of
the reasons why we used the ARI is that it can show the performance of the
clustering qualitatively since we already have the original information of the
network However, with observational data we do not know the true networks.
Consequence, we can only track the number of detected networks changing with
both linking lengths and with the minimum number of member galaxies from 2 to
5, and comment on the results.

For all four clustering algorithms, the HR4 simulation results match well with
the observations within cosmic variance, especially for $n\geq 5$ at
$\ell\ga\langle d_{\rm particle}\rangle=0.5\mathinner{h^{-1}{\rm Mpc}}$. On
the other hand, HR4 tends to underestimate the number of networks for a
smaller minimum number of member galaxies and/or smaller linking length
$\ell\la 0.5\mathinner{h^{-1}{\rm Mpc}}$. This may mean that, despite the
agreement with the observation in terms of 2pCF below $1\mathinner{h^{-1}{\rm
Mpc}}$-scale, some disagreements exist between HR4 and observation in terms of
the higher-order statistics in smaller scales than the particle mean
separation scale.

One notable feature of the MGS algorithm is that it does not create a single
giant network for large linking lengths. Instead, the algorithm identifies a
number of smaller networks, even for large linking lengths. This is in
contrast to the other three algorithms, which all create a single giant
network for large linking lengths. This difference highlights the ability of
the MGS algorithm to accurately identify networks in data with a complex
distribution of galaxies.

Figure 11: The number of member galaxies for the 1st, 2nd, 3rd and 4th richest
networks in D4-SDSS (thick lines) and in D4-HR4 as a function of linking
length. For D4-HR4, the average values and the ranges between minimums and the
maximums of 4 data samples are drawn as thin lines and error bars. Color
indicates the different clustering algorithms.

Fig. 11 shows the number of member galaxies for the 1st to 4th richest
networks identified by each algorithm in D4-SDSS and D4-HR4. Similar to Fig.
10, both results from the simulation and observation data match well with each
other. The top left panel of the figure shows the shape of the richest network
for each algorithm. As the linking length increases over certain value, the
richest network identified by the FoF, MST, and DBSCAN algorithms contains all
of the galaxies in the data, while the MGS algorithm identifies a network with
only a portion of the galaxies. This indicates that the MGS algorithm is able
to identify multiple networks even for large linking lengths, while the other
algorithms merge all of the galaxies into a single giant network.

The main difference between the MGS algorithm and the other three algorithms
becomes particularly clear when examining the number of member galaxies in the
2nd to 4th richest networks (upper right and bottom panels of Fig. 11). As the
linking length increases, the number of member galaxies in these networks
identified by the FoF, MST, and DBSCAN algorithms decreases to zero. This is
because the first richest network identified by these algorithms took all the
galaxies in the data, leaving no galaxies to be considered for further
clustering. In contrast, the MGS algorithm is able to identify multiple
networks even with large linking lengths as the richest network does not
monopolize all galaxies. This demonstrates the ability of the MGS algorithm to
accurately identify networks in data with a complex distribution of galaxies.

Figure 12: Top 30 richest networks (colors) found by the MGS in the D4-SDSS
galaxies (gray dots). The observer is located at the origin. The linking
length is $\ell_{\rm MGS}=10\mathinner{h^{-1}{\rm Mpc}}$, where all D4-SDSS
galaxies fall into a single giant network in all other three benchmark
algorithms (see Fig. 11). Note that, even in such a large linking length, none
of the 30 richest networks suffers percolation.

Fig. 12 shows the 30 richest networks found by the MGS algorithm in the
D4-SDSS data with a linking length $\ell_{\rm MGS}=10\mathinner{h^{-1}{\rm
Mpc}}$. While such a large choice of linking length makes a single giant
network in all other three benchmark algorithms (see Fig. 11), none of the 30
networks suffer percolation. All 30 richest networks are well-separated and
have some even distribution of galaxies in the XY-plane (that is, the
tangential plane). On the other hand, most of the networks have somewhat
elongated features in the line-of-sight direction, which clearly shows the
Finger-of-God effect due to the RSD. Therefore, although this needs further
inspection, we consider that the 30 richest networks found by the MGS
algorithm could be the actual large structures similar to galaxy
(super)clusters in real space. We stress that, at this stage, the MGS
algorithm is not necessarily finding gravitationally bound clusters, but
rather simply looking at distinct over-densities in the galaxy distribution.
In future work we will explore the physical properties of these networks and
evaluate their usefulness in cosmological analyses.

##  5 Conclusions

The MulGuisin (MGS) algorithm is a powerful technique for identifying networks
in data from astrophysical simulations and observations. It consistently
produces results closer to those inferred from human visual inspection. In
comparison to other clustering algorithms, such as the friends-of-friends
(FoF) algorithm, the minimum spanning tree (MST) algorithm, and the DBSCAN
algorithm, the MGS algorithm has several advantages. The MGS algorithm is able
to take into consideration the local density and is able to accurately
identify networks even in complex data sets with a large number of galaxies.
In contrast, the FoF, MST, and DBSCAN algorithms often merge networks into a
single giant network for large linking lengths, losing the ability to
accurately identify individual networks. This characteristic of the MGS
algorithm is particularly important for analyzing data from astrophysical
simulations and observations.

In this proof of concept work we have shown that the jet-finding algorithm MGS
can be applied to mock Galaxy data resulting in reliable network
identification. However the identification of networks in real observation is
a difficult issue due to survey incompleteness, selection effects, redshift-
space distortions, etc. In future work we will test MGS in the presence of
realistic observational systematic effects.

MGS also provides auxiliary topological information such as the number and
length of connections for each galaxy. In future work we will explore the use
of this enhanced information in testing or constraining cosmological models.

## Acknowledgements

The authors thank Changbom Park, Dongsu Bak, and Ena Choi for helpful
discussions.

This research was supported by Basic Science Research Program through the
National Research Foundation of Korea(NRF) funded by the Ministry of Education
(2018R1A6A1A06024977) for Young Ju. C.G.S is support via the Basic Science
Research Program from the National Research Foundation of South Korea (NRF)
funded by the Ministry of Education (2020R1I1A1A01073494). S.E.H. was
supported by the project 우주거대구조를 이용한 암흑우주 연구 (“Understanding Dark Universe
Using Large Scale Structure of the Universe”), funded by the Ministry of
Science. This work was supported by the 2023 Research Fund of the University
of Seoul for Inkyu Park.

This work was supported by the Supercomputing Center/Korea Institute of
Science and Technology Information, with supercomputing resources including
technical support (KSC-2013-G2-003), and the simulation data were transferred
through a high-speed network provided by KREONET/GLORIAD.

We would like to acknowledge the Korean Astronomy and Machine Learning (KAML)
meeting series for providing a forum for fruitful discussions on the
intersection of astronomy and machine learning.

Funding for the SDSS and SDSS-II has been provided by the Alfred P. Sloan
Foundation, the Participating Institutions, the National Science Foundation,
the US Department of Energy, the National Aeronautics and Space
Administration, the Japanese Monbukagakusho, the Max Planck Society, and the
Higher Education Funding Council for England. The SDSS website is
http://www.sdss.org/.

The SDSS is managed by the Astrophysical Research Consortium for the
Participating Institutions. The Participating Institutions are the American
Museum of Natural History, Astrophysical Institute Potsdam, University of
Basel, University of Cambridge, Case Western Reserve University, University of
Chicago, Drexel University, Fermilab, the Institute for Advanced Study, the
Japan Participation Group, Johns Hopkins University, the Joint Institute for
Nuclear Astrophysics, the Kavli Institute for Particle Astrophysics and
Cosmology, the Korean Scientist Group, the Chinese Academy of Sciences
(LAMOST), Los Alamos National Laboratory, Max Planck Institute for Astronomy
(MPIA), the Max Planck Institute for Astrophysics (MPA), New Mexico State
University, Ohio State University, University of Pittsburgh, University of
Portsmouth, Princeton University, the US Naval Observatory, and the University
of Washington.

## Data Availability

The up-to-date MulGuisin algorithm can be downloaded at
https://github.com/youngju20/Mulguisin.

## References

  * Abazajian et al. (2009) Abazajian K. N., et al., 2009, ApJS, 182, 543
  * Alcock & Paczynski (1979) Alcock C., Paczynski B., 1979, Nature, 281, 358
  * Appleby et al. (2017) Appleby S., Park C., Hong S. E., Kim J., 2017, ApJ, 836, 45
  * Appleby et al. (2018) Appleby S., Chingangbam P., Park C., Hong S. E., Kim J., Ganesan V., 2018, ApJ, 858, 87
  * Bae & Woo (2014) Bae H.-J., Woo J.-H., 2014, ApJ, 795, 30
  * Ballinger et al. (1996) Ballinger W. E., Peacock J. A., Heavens A. F., 1996, MNRAS, 282, 877
  * Barrow et al. (1985) Barrow J. D., Bhavsar S. P., Sonoda D. H., 1985, Monthly Notices of the Royal Astronomical Society, 216, 17
  * Berlind & Weinberg (2002) Berlind A. A., Weinberg D. H., 2002, The Astrophysical Journal, 575, 587
  * Blanton et al. (2005) Blanton M. R., et al., 2005, AJ, 129, 2562
  * Blumenthal et al. (1984) Blumenthal G. R., Faber S. M., Primack J. R., Rees M. J., 1984, Nature, 311, 517
  * Borůvka (1926) Borůvka O., 1926, Plàce mor. přírodově d. spol. v Brně III (3), pp 37–58
  * Bosman et al. (1998) Bosman M., Park I., Corbal M., Costanzo D., Lami S., Paoletti R., Azuelos G., Strahl K., 1998, ATLAS Software, 98, 038
  * Choi et al. (2010a) Choi Y.-Y., Han D.-H., Kim S. S., 2010a, Journal of Korean Astronomical Society, 43, 191
  * Choi et al. (2010b) Choi Y.-Y., Park C., Kim J., Gott J. Richard I., Weinberg D. H., Vogeley M. S., Kim S. S., SDSS Collaboration 2010b, ApJS, 190, 181
  * Dall & Christensen (2002) Dall J., Christensen M., 2002, Phys. Rev. E, 66, 016121
  * Davis et al. (1985) Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, ApJ, 292, 371
  * DeRose et al. (2019) DeRose J., et al., 2019, ApJ, 875, 69
  * Diemer (2018) Diemer B., 2018, The Astrophysical Journal Supplement Series, 239, 35
  * Duarte & Mamon (2014) Duarte M., Mamon G. A., 2014, Monthly Notices of the Royal Astronomical Society, 440, 1763
  * Dunkley et al. (2009) Dunkley J., et al., 2009, ApJS, 180, 306
  * Ester et al. (1996) Ester M., Kriegel H.-P., Sander J., Xu X., et al., 1996, in kdd. pp 226–231
  * Fall & Efstathiou (1980) Fall S. M., Efstathiou G., 1980, MNRAS, 193, 189
  * Farrens et al. (2011) Farrens S., Abdalla F. B., Cypriano E. S., Sabiu C., Blake C., 2011, Mon. Not. Roy. Astron. Soc., 417, 1402
  * Gifford et al. (2013) Gifford D., Miller C., Kern N., 2013, ApJ, 773, 116
  * Gott et al. (1986) Gott J. Richard I., Melott A. L., Dickinson M., 1986, ApJ, 306, 341
  * Guo et al. (2012) Guo H., Zehavi I., Zheng Z., 2012, ApJ, 756, 127
  * Hearin et al. (2017) Hearin A. P., et al., 2017, The Astronomical Journal, 154, 190
  * Hong et al. (2016) Hong S. E., Park C., Kim J., 2016, ApJ, 823, 103
  * Hoyle et al. (2012) Hoyle F., Vogeley M. S., Pan D., 2012, MNRAS, 426, 3041
  * Hubert & Arabie (1985) Hubert L., Arabie P., 1985, Google Scholar, pp 193–128
  * Huchra & Geller (1982) Huchra J. P., Geller M. J., 1982, The Astrophysical Journal, 257, 423
  * Hwang et al. (2012) Hwang H. S., Park C., Elbaz D., Choi Y. Y., 2012, A&A, 538, A15
  * Jiang et al. (2008) Jiang C. Y., Jing Y. P., Faltenbacher A., Lin W. P., Li C., 2008, ApJ, 675, 1095
  * Kim et al. (2015) Kim J., Park C., L’Huillier B., Hong S. E., 2015, Journal of Korean Astronomical Society, 48, 213
  * Knebe et al. (2011) Knebe A., et al., 2011, MNRAS, 415, 2293
  * Kravtsov et al. (2004) Kravtsov A. V., Berlind A. A., Wechsler R. H., Klypin A. A., Gottlober S., Allgood B., Primack J. R., 2004, The Astrophysical Journal, 609, 35
  * Krzewina & Saslaw (1996) Krzewina L. G., Saslaw W. C., 1996, Monthly Notices of the Royal Astronomical Society, 278, 869
  * Lee et al. (2012a) Lee G.-H., Park C., Lee M. G., Choi Y.-Y., 2012a, ApJ, 745, 125
  * Lee et al. (2012b) Lee G.-H., Woo J.-H., Lee M. G., Hwang H. S., Lee J. C., Sohn J., Lee J. H., 2012b, ApJ, 750, 141
  * Li et al. (2014) Li X.-D., Park C., Forero-Romero J. E., Kim J., 2014, ApJ, 796, 137
  * Limousin et al. (2013) Limousin M., Morandi A., Sereno M., Meneghetti M., Ettori S., Bartelmann M., Verdugo T., 2013, Space Sci. Rev., 177, 155
  * More et al. (2011) More S., Kravtsov A. V., Dalal N., Gottlöber S., 2011, ApJS, 195, 4
  * Morey & Agresti (1984) Morey L. C., Agresti A., 1984, Educational and Psychological Measurement, 44, 33
  * Naidoo (2019) Naidoo K., 2019, Journal of Open Source Software, 4, 1721
  * Naidoo et al. (2020) Naidoo K., Whiteway L., Massara E., Gualdi D., Lahav O., Viel M., Gil-Marún H., Font-Ribera A., 2020, Monthly Notices of the Royal Astronomical Society, 491, 1709
  * Navarro et al. (1996) Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
  * Pan et al. (2012) Pan D. C., Vogeley M. S., Hoyle F., Choi Y.-Y., Park C., 2012, MNRAS, 421, 926
  * Park & Gott (1991) Park C., Gott J. R. I., 1991, ApJ, 378, 457
  * Park & Kim (2010) Park C., Kim Y.-R., 2010, ApJ, 715, L185
  * Park et al. (2012) Park C., Choi Y.-Y., Kim J., Gott J. Richard I., Kim S. S., Kim K.-S., 2012, ApJ, 759, L7
  * Park et al. (2019) Park H., Park C., Sabiu C. G., Li X.-d., Hong S. E., Kim J., Tonegawa M., Zheng Y., 2019, ApJ, 881, 146
  * Pedregosa et al. (2011) Pedregosa F., et al., 2011, Journal of Machine Learning Research, 12, 2825
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A13
  * Press & Schechter (1974) Press W. H., Schechter P., 1974, ApJ, 187, 425
  * Rand (1971) Rand W. M., 1971, Journal of the American Statistical association, 66, 846
  * Reid et al. (2014) Reid B. A., Seo H.-J., Leauthaud A., Tinker J. L., White M., 2014, MNRAS, 444, 476
  * Rykoff et al. (2014) Rykoff E. S., et al., 2014, ApJ, 785, 104
  * Sander et al. (2017) Sander J., Ester M., et al., 2017, ACM Transactions on Database Systems, 42, 1
  * Serra & Diaferio (2013) Serra A. L., Diaferio A., 2013, ApJ, 768, 116
  * Sheth (1996) Sheth R. K., 1996, MNRAS, 279, 1310
  * Tago et al. (2008) Tago E., Einasto J., Saar E., Tempel E., Einasto M., Vennik J., Müller V., 2008, Astronomy & Astrophysics, 479, 927
  * Tempel et al. (2016) Tempel E., Kipper R., Tamm A., Gramann M., Einasto M., Sepp T., Tuvikene T., 2016, Astronomy & Astrophysics, 588, A14
  * Tonegawa et al. (2020) Tonegawa M., Park C., Zheng Y., Park H., Hong S. E., Hwang H. S., Kim J., 2020, ApJ, 897, 17
  * White & Rees (1978) White S. D. M., Rees M. J., 1978, MNRAS, 183, 341
  * Yang et al. (2007) Yang X., Mo H. J., van den Bosch F. C., Pasquali A., Li C., Barden M., 2007, ApJ, 671, 153
  * Zehavi et al. (2002) Zehavi I., et al., 2002, ApJ, 571, 172
  * Zehavi et al. (2011) Zehavi I., et al., 2011, ApJ, 736, 59