Unsupervised Machine Learning in Astronomy

1. Introduction

Machine learning (ML) has been applied to various analyses of astronomical data, such as analyzing spectral data e.g., [1,2,3], catalogs e.g., [4,5,6], light curves e.g., [7,8,9], and images e.g., [10,11,12]. ML, a subfield of artificial intelligence, aims to mimic the human brain using computers. Rather than manually coding every step, researchers use existing ML models to analyze data. In simple terms, ML performs tasks with general guidance rather than detailed instructions [13]. Unlike traditional programming, ML algorithms often involve iterative steps that are not easily described by equations.

The growing significance of ML arises from the rapidly increasing volume and the complexity of astronomical data collected using progressively advanced instruments [14]. With the advanced instruments, innumerable high-resolution, high-dimensional data sets are collected. ML offers efficient and objective solutions for analyzing such large datasets. Therefore, ML has become increasingly popular among astronomers.

There are two types of ML techniques: unsupervised and supervised. Unsupervised ML conducts exploratory data analysis, discovering unknown data features [15], without any prior information on classification. In comparison, supervised ML requires a labeled dataset (i.e., a training set), where the features are known [14]. Supervised ML learns from this human-labeled dataset and predictions on new data within the same properties [14,16]. The project focuses on unsupervised ML algorithms, which are sometimes considered more helpful for scientific research as they are not limited by present knowledge and can be used to extract new knowledge [14].

In the review, we classify the unsupervised ML algorithms into three categories: clustering, dimensionality reduction, and neural network. Clustering and dimensionality reduction are often considered objectives, whereas neural network is considered a model that extracts information from the data, and the information can be used to perform tasks such as clustering and dimensionality reduction.

Clustering refers to finding the concentration of multivariable data points [15]. In simpler words, clustering groups objects so that those in the same group are more similar to each other than to those in different groups.

Dimensionality reduction selects or constructs a subset of features that best describe the data, reducing the number of features [14]. It retains essential information while discarding trivial information [15]. An important branch of dimensionality reduction - manifold learning - performs non-linear reduction to unfold the surface of data and reveal its underlying structure [15].

A neural network model is designed to mimic the structure and function of the human brain [15]. A neural network consists of multiple interconnected neurons or layers of neurons, where each neuron receives inputs and transmits outputs. A shallow neural network has one or two hidden layers, while deep learning models have three or more. While deep learning can be used for supervised ML, this review focuses on unsupervised neural networks.

This review aims to provide an overview of unsupervised ML algorithms used in astronomical research to analyze data, classified by the type of tasks. The review may serve as an up-to-date, comprehensive manual on the topic, tailored for astronomy researchers new to ML.

2. Dimensionality Reduction

This section introduces various dimensionality reduction algorithms, including principal component analysis, multi-dimensional scaling and isometric feature mapping, locally linear embedding, and t-distributed stochastic neighbor embedding. These algorithms project high-dimensional data into lower dimensions by identifying linear or nonlinear structures that preserve essential information and, in some cases, by studying the underlying manifold of the data.

2.1. Principal Component Analysis (PCA) and Kernal PCA

Principal component analysis (PCA), as suggested by its name, is a dimensionality reduction method that focuses on the principal components of the multivariate data set. PCA was invented by Pearson [17] and developed by Hotelling [18,19]. PCA performs singular value decomposition (SVD) and can be seen as a rotation and projection of the data set that maximizes the variance, making the data features more significant [15]. PCA constructs a covariance matrix of the data set, and the orthonormal eigenvectors are the principal components (i.e., the axes) [14]. The principal components are identified sequentially: the first maximizes the variance, the second is orthogonal to the first while maximizing the residual variance, and the subsequent components are orthogonal to all prior components [15]. The principal components are linearly uncorrelated, so applying PCA removes the correlation between multiple dimensions, simplifying the data. The first few components convey most of the information [16]. Therefore, when all principal components are used, the full information of the data set is preserved; when k principal components are used, the data are reduced to k-dimension. Figure 1 is a demonstration of dimensionality reduction with PCA, showing how the two principal components are selected given the data set. Figure 2 illustrates an example of PCA applied to images, showing how the reconstructed images retain the main features of the original data. The number n on the left indicates the number of principal components used for the reconstruction. A smaller number of components corresponds to greater compression in feature space, which results in a blurrier reconstructed image.

There are some shortcomings in PCA. Firstly, pre-processing, i.e., treatment of the data set before applying the algorithm, is required to produce more informative results [15]. PCA is sensitive to outliers, so outliers need to be removed. Then, due to the property of PCA, the data needs to be normalized beforehand, similar to K-means (see §3.2) and some other algorithms. For instance, z-score normalization, sometimes referred to as feature scaling, can be applied. The equation is $z=(x-\mu )/\sigma$, where z is the z-score, x is the value being evaluated, and $\mu$ and $\sigma$ are respectively the mean and standard deviation of the data. Secondly, PCA is a linear decomposition of data, thus not applicable in some cases (e.g., when effects are multiplicative) [14]. When applied to analyze a nonlinear data set, PCA may fail. Thirdly, exact PCA only supports batch processing, requiring all data to fit in the main memory. Therefore, incremental PCA is developed to support minibatch processing but is relatively less accurate.

As mentioned, PCA is a linear method. Therefore, to analyze data sets that are not linearly separable, kernel PCA [22] is introduced. Compared to PCA, kernel PCA is able to make a non-linear projection of the data points, thus providing a clearer presentation of information by unfolding the dataset. When mapping the kernelized PCA back to the original feature space, there will be some small differences even if the number of components defined is the same as the number of original features. The difference can be reduced by imposing a different built-in function or revising the code, as suggested by Pedregosa et al. [23].

PCA has been applied to dimensionality reduction of spectral data e.g., [24,25,26,27,28], light curves e.g., [29], catalogs e.g., [30], and images [11].

2.2. Multi-dimensional Scaling and Isometric Feature Mapping

Multi-dimensional scaling (MDS) [31] is a dimensionality reduction algorithm frequently compared to PCA. James O. Ramsay [32] gives the core theory behind MDS. MDS aims to preserve the disparities, which are the pairwise distances computed between data points. When given a dataset, MDS first finds the pairwise distances between all pairs of two points, then reconstructs a low-dimensional dataset that minimizes stress (a.k.a. the error). The stress is the summation of the differences between the original distance and the distances in the lower-dimensional space. There are various ways to define the disparities. For example, the MDS metric, a.k.a. the absolute MDS, defines disparity as a factor of the absolute distance between two points. In contrast, the non-metric MDS [33,34] forces the rank order of the distances between pairs to be the same. That is, if one pair of points is farther apart compared to another pair of points, the same relationship remains in the embedding space.

Isometric feature mapping (Isomap) [35], an extension of MDS, is also used for nonlinear dimensionality reduction. Instead of using Euclidean distance as in PCA and MDS, Isomap uses geodesic distance. Based on the Euclidean distance, Isomap connects each point to its nearest neighbors, constructing a neighborhood graph. The geodesic distance is then approximated as the shortest path between two points along this graph. In other words, Isomap approximates the geodesic curves lying within the manifold and computes the distance along the geodesic curves [15]. While MDS considers the distances between all pairs of points to define the shape, Isomap only uses the distances between neighboring points and sets the distances between any other two points to be 0, thus unfolding the manifold [15]. Then, Isomap applies MDS to this new shape. Figure 3 (a) shows an example of how a three-dimensional ‘Swiss roll’ manifold is flattened to two dimensions by Isomap.

Bu et al. [36] points out that Isomap is more efficient than PCA in feature extraction for spectral classification. However, Ivezić et al. [15] states that Isomap is more computationally expensive compared to locally linear embedding (see §2.3). The book also lists various algorithms that can reduce the amount of computation, such as the Floyd–Warshall algorithm [37] and the Dijkstra algorithm [38].

MDS has been used to reduce the dimension of catalogs e.g., [12], while Isomap has been used to reduce spectral data e.g., [36,39,40] and latent space generated by neural network (see §4.2) e.g., [41].

2.3. Locally Linear Embedding

Locally linear embedding LLE [42] is an algorithm for dimensionality reduction that preserves the data’s local geometry [15]. For each data point, LLE identifies its k nearest neighbors and produces a set of weights that can be applied to the neighbors to best reconstruct the data point. The weight indicates the geometry defined by the data point and its nearest neighbors. The weight matrix W is find by minimizing the error ${\mathcal{E}}_1\left(W\right)={|X-WX|}^2$, where X denotes the original data set [15]. From the weight matrix, LLE finds the low-dimensional embedding Y by minimizing the error ${\mathcal{E}}_2\left(Y\right)={|Y-WY|}^2$ [15]. The solution to the two equations can be find by imposing efficient linear algebra techniques: some computations are made to find another matrix $C_W\equiv {(I-W)}^T(I-W)$, and eigenvalue decomposition is performed on $C_W$, producing the low-dimensional embedding. Figure 3 (b) shows an example of how a three-dimensional ‘Swiss roll’ manifold is flattened to two dimensions by LLE.

A shortcoming of LLE is that the direct eigenvalue decomposition becomes expensive for large data [15]. According to Ivezić et al. [15], this problem can be circumvented by using iterative methods, such as Arnoldi decomposition, which is available in the Fortran package ARPACK [43]. LLE is sometimes compared with PCA. On the one hand, LLE does not allow the projection of new data, which would impact the weight matrix and subsequent computations. This means that when LLE is applied to reduce new data of the same type, the training process needs to be repeated, unlike some other methods that allow directly applying a trained template to new data. On the other hand, Vanderplas and Connolly [44] shows that LLE leads to improved classification, although Bu et al. [45] points out that the usage of LLE is more specific than PCA, thus more limited.

LLE has been applied to dimensionality reduction of spectral data e.g., [44,46,47,48] and light curves e.g., [7,49].

2.4. t-distributed Stochastic Neighbor Embedding

t-distributed Stochastic Neighbor Embedding (t-SNE) is another widely used algorithm for dimensionality reduction. Stochastic Neighbor Embedding is the basis of t-SNE, developed by Hinton and Roweis [50]. The t-distribution variant is added by van der Maaten and Hinton [51]. t-SNE converts the affinities - similarities between two points, usually measured by the pairwise distances in a scattered space - of the data points into Gaussian joint probabilities [23]. Figure 4 illustrates how t-SNE calculates affinity, using a three-dimensional data set shown in (a). For instance, the distance from any dot to the darkest black dot is computed, and the t-SNE finds the conditional Gaussian probabilities, as shown in (b). Perplexity, a hyperparameter that reflects the effective number of nearest neighbors considered when computing conditional probabilities, determines the width of the Gaussian distribution around each point. Because of the normalization of the Gaussian distribution, the probabilities between any points a and b can differ depending on whether the perspective is from a to b or from b to a. The joint probability is typically defined as the average of these two probabilities. The probabilities are further represented by Student’s t-distributions [23]. When high-dimensional data is projected into a low-dimensional space, the probability distribution is maximally preserved [14]. That is, t-SNE uses joint Gaussian distribution to model the likelihood of data in the high-dimensional space, and the Gaussian distribution is mapped to a Student’s t-distribution in the low-dimensional space [16]. Figure 5 illustrates an example of t-SNE applied to spectral data. Each spectrum can be represented as a vector whose dimensionality equals the total number of wavelength bins. The color of each dot indicates the cluster it belongs to, as determined by a given clustering algorithm, while the pink dots represent intermediate data points. The t-SNE projection reveals that the circled groups of pink data points lie near the boundary between the two main clusters, suggesting the presence of possible intermediate-type subgroups.

As discussed, Isomap and LLE can learn a single continuous manifold. However, in many cases, there is more than one manifold. Due to its design, t-SNE can learn different manifolds in the dataset. t-SNE also addresses the issue of overcrowding at the center - a common problem in many dimensionality reduction algorithms, where points become densely packed in the center - by using a t-distribution, which has heavier tails than a Gaussian distribution [16]. This allows distant points to be spread out more effectively in the lower-dimensional space. Nonetheless, t-SNE has some disadvantages. Similar to LLE, t-SNE does not allow the projection of new data points once the manifold has been learned. t-SNE is non-deterministic, which means it generates different results each time it runs. t-SNE is also very computationally expensive. To accelerate computation, the Barnes-Hut approximation [53] may be adopted, but the embedding manifolds applicable are limited to two or three dimensions. t-SNE may not preserve the global structure of the data set. To solve this problem, one may select initial points by PCA.

t-SNE has been applied to reduce the dimensionality of catalogs e.g., which are respectively quasars parameters, line ratios, and chemical abundance [54,55,56], spectral data [57], photometry data [58], and light curves [8,59].

2.5. Examples of Applications of Different Dimensionality Reduction Algorithms

In this section, we present examples of applying different dimensionality reduction algorithms to real observational datasets in astronomy. The data set consists of five-dimensional astrometry data for 1,254 stars, including right ascension, declination, distance, proper motion in right ascension, and proper motion in declination. The data is provided by SIMBAD [60]. We used Strasbourg Astronomical Data Center CDS, [61] for the criteria query of data. The criteria we used was $55°<RA<70°$, $20°<DEC<35°$, $100\mathrm{pc}<\mathrm{Distance}<170\mathrm{pc}$, and Spectral types being dimmer or equal to ‘O’. The data is pre-processed by normalizing, removing outliers, and then re-normalizing to adjust the scale based solely on the inlier data, retaining 1,113 data points.

Figure 6 shows the results of the dimensionality reduction algorithms, including PCA, LLE, non-metric MDS, and t-SNE. The five-dimensional data set is clustered with DBSCAN (see §3.4) before being reduced to two dimensions. PCA is linear, while the other three algorithms are non-linear.

3. Clustering

Given a data set, we may be interested in classifying it into groups so that each group contains similar data. Unsupervised clustering algorithms are designed to accomplish this task without the need for human-provided labels. This section introduces various clustering algorithms, including the Gaussian mixture model, the K-means algorithm, hierarchical clustering, density-based spatial clustering of applications with noise and its hierarchical variant, and fuzzy C-means clustering. These algorithms group data sets into clusters by either iteratively estimating functions that represent the clusters, expanding clusters outward from core data points, or progressively merging or splitting clusters.

3.1. Gaussian Mixture Model

The Gaussian mixture model (GMM) is widely applied to group objects or data into clusters. GMM is a probability model that assumes the data distribution follows the weighted summation of multiple Gaussian functions, where each function has a weight, and the sum of the weights is 1. Each Gaussian function can be seen as a component. Therefore, it is possible to cluster the objects based on which weighted Gaussian function has the highest contribution at that point. The expectation-maximization (EM) algorithm is the iterative process of fitting Gaussian functions to the data. The EM algorithm first computes the probability of generating each data point using Gaussian functions with random parameters. Then, the EM algorithm changes the parameters to maximize the probability. GMM assumes the clusters to be convex-shaped, i.e., each cluster has a single center and follows a relatively ellipsoidal distribution around it. Figure 7 shows an example in which two components are considered for clustering. The lines can be seen as the contour line of the sum of the two Gaussian functions.

GMM assumes that the number of components is known. However, in most applications, the number of components is unknown. As a result, the users may either apply the Bayesian information criterion BIC, [62] to determine the number of components or use the variational Bayesian Gaussian mixture model (VBGMM), which does not fix the number of components. The BIC is a score computed from a grid search over different numbers of components and the shapes of their distributions (i.e., the types of covariance). The combination with the lowest BIC score indicates the best fit to the data distribution and is used for GMM. The VBGMM requires more hyperparameters than EM. There are different types of weight concentration prior, which is an important parameter for VBGMM. The VBGMM using a Dirichlet distribution prior fixes the maximum number of components and uses a concentration parameter that controls the weighting of the components: the lower the concentration parameter, the more weight is placed on fewer components [23]. On the other hand, the VBGMM using a Dirichlet process prior (i.e., DPGMM) may have infinite components, and the concentration parameter is used to restrain the number of components likewise [23].

GMM has been applied to cluster photometric data e.g., [63,64], spectroscopic data e.g., [65], catalogs e.g., [66,67], and scattered data e.g., [68].

3.2. K-means

The K-means algorithm was invented by MacQueen [69] to partition N-dimensional objects into k clusters. The mechanism of the K-means algorithm is explained as follows. Firstly, centroid initialization is applied. From the data, k initial centroids are picked, and each elements are assigned to be in the same cluster as the closest centroid. Then, a new centroid for each cluster is computed. Given the position of the new centroids, repeat the assignment and recompute the new centroids. K-means will eventually converge, leaving centroids that produce minimal in-cluster sum-of-squares (i.e., the sum of the squared distances between the centroid and elements in that cluster).

The K-means algorithm has some disadvantages. Similar to the Gaussian mixture model (§3.1), K-means assumes prior knowledge of the number of clusters. There are ways to avoid the problem, such as trying different numbers of clusters and selecting the one where the distortion curve begins to level off, indicating diminishing gain from adding more clusters [6], or using a large number of clusters and discarding the small clusters if the data permits [1]. Another disadvantage of the K-means algorithm is its high dependence on the centroid initialization. One approach to solve this problem is to repeat the computation for different initializations and use the case that produces the smallest in-cluster sum-of-squares, as demonstrated in [3]. Another approach is to use a variant of the K-means algorithm, called the K-means++ algorithm [70], which forces the initial guess of the centroids to be distant from each other to improve the clustering. K-means also assumes the clusters to be convex-shaped.

The K-means algorithm has been applied to various types of data for clustering, including scattered data from catalogs e.g., [6], spectral data e.g., [1,2,3], and polarimetric data [71].

3.3. Hierarchical Clustering (HC)

Hierarchical clustering (HC) is a clustering algorithm that identifies clusters across all scales without requiring a specified number of clusters [15]. HC can be performed in a top-down (a.k.a. divisive) or bottom-up a.k.a. agglomerative; [72] approach. A diagram demonstrating the bottom-up clustering dendrogram is Figure 8. In the bottom-up approach, the N elements start as N independent clusters, each containing a single element. The two nearest clusters are then merged, reducing the number of clusters from N to $N-1$. The distance between clusters is calculated in various ways, generating significantly different results. One example is to sum the distances between all possible pairs of points, with one point from each cluster, and divide by the product of the two numbers of elements from the two clusters [15]. The merging procedure is repeated until there is only one cluster, containing all N elements. The top-down approach is the reverse of the bottom-up approach. Instead of merging clusters, this method divides a single cluster into two at each step. By imposing HC, all possible clustering with all possible numbers of clusters are generated, and the user can choose the level (i.e., the number of clusters) for investigation.

An advantage of HC is that the computation does not need to be repeated if different numbers of clusters are considered. Another advantage is that HC does not assume the shape of clusters to be convex, unlike GMM and K-means algorithms.

HC has been applied to cluster bivariate data e.g., chemical abundances and positions of stars in [73,74], higher-dimensional scattered data e.g., [75], light curves e.g., [9], and spectral data e.g., [76].

3.4. Density-Based Spatial Clustering of Applications with Noise

Density-based spatial clustering of applications with noise DBSCAN, [77] is another clustering algorithm that discovers clusters of all shapes and does not require a specified number of clusters. DBSCAN divides the areas into high-density areas and low-density areas. High-density areas are where the cluster lies. Core samples are picked from high-density areas. We define the core sample as a sample that has a minimum of k samples at most distance s from itself. Samples within a distance s from a core sample are considered neighbors of the sample. Then, the same requirement is applied to find the core samples among the neighbors. The samples are considered to be in the same cluster as their neighboring core sample. Then, for each additional core sample, we find its neighboring core samples. The steps are performed recursively to cluster the samples.

DBSCAN has some advantages and disadvantages. As discussed above, an advantage of DBSCAN is that it does not assume the clusters to have convex shapes, thus discovering clusters with arbitrary shapes. Another advantage is that DBSCAN can filter out the outliers, which are non-core samples that do not neighbor any core sample. On the other hand, one crucial disadvantage of DBSCAN is that the results are greatly impacted by the parameters k and s, especially s. For high-density data, DBSCAN may require a higher k for better clustering. s is highly dependent on the data: setting s too small would lead to the fringes of clusters being recognized as outliers, while setting s too large could merge clusters. However, there are various ways to determine k and s, e.g., [78,79,80,81]. Another disadvantage of DBSCAN is the single density threshold conveyed by the fixed k and s, which means DBSCAN may not be useful when the clusters have different densities. The disadvantage could be avoided using hierarchical DBSCAN, which is introduced in §3.5.

DBSCAN has been applied to five-dimensional scattered data from catalogs e.g., position and motion taken from GAIA catalogs by [78,82,83], three-dimensional positional scattered data e.g., positions of stars, taken from GAIA catalogs, by [79], spectral data [80,81], and images [10].

3.5. Hierarchical Density-Based Spatial Clustering of Applications with Noise

Hierarchical density-based spatial clustering of applications with noise HDBSCAN, [84] is the hierarchical extension of DBSCAN, as indicated by the name. HDBSCAN is also used for clustering, especially for globally inhomogeneous data, whereas DBSCAN assumes the data to be homogeneous. As discussed above, DBSCAN imposes a single density threshold to define clusters. As a result, DBSCAN may not be applicable when the clusters have different densities. HDBSCAN solves this problem by fixing the minimum number of samples k and considering all possible distance s (for more information on k and s, refer to §3.4). First, HDBSCAN defines the core distance of a point as the distance to its nearest k-th point. For all pairs of points p and q, HDBSCAN defines the mutual reachability distance as the maximum of the core distances of the two points and the distance between them and thus transforms the graph so that each pair of data points is separated by their mutual reachability distance. Then, HDBSCAN finds the minimum spanning tree in the new graph. From the minimum spanning tree, HDBSCAN uses HC to find all possible clusterings.

HDBSCAN has all the advantages of DBSCAN, but it also has two additional advantages: HDBSCAN does not apply the same density threshold to all clusters, and the computation does not need to be repeated to consider different numbers of clusters. HDBSCAN eliminates the use of s and instead has a new parameter - the minimum size of a cluster. In some cases, this parameter may be easier to set than s, since it is basically asking what size of a group of data you would consider a cluster [85].

HDBSCAN has been applied to spectroscopic and photometric data e.g., [86,87,88], light curves e.g., [89], astrometric data e.g., [90], and other catalogs e.g., [91,92].

3.6. Fuzzy C-means Clustering

Fuzzy C-means clustering (FCC, a.k.a. C-means clustering), being the most widely applied fuzzy clustering (a.k.a. soft clustering) algorithm, was invented by Dunn [93] and further developed by Bezdek [94]. Soft clustering means that the algorithms may not assign a data point to a single cluster. Instead, the point may be assigned to multiple clusters with corresponding membership grades between 0 and 1. That is, a point at the edge of a cluster may have a smaller membership grade for that cluster (e.g., 0.1) compared to a point at the center of the cluster (e.g., 0.95). Turning to FCC, it is very similar to the K-means algorithm mentioned in §3.2. The only difference is that K-means clustering imposes hard clustering, while C-means clustering imposes soft clustering. In fact, K-means clustering is sometimes referred to as hard C-means clustering [95]. Figure 9 shows an application of FCC to two-dimensional scattered data from the Sloan Moving Object Catalog.

The advantages and disadvantages are very similar to those of K-means clustering. One shared disadvantage is the need for prior knowledge about the number of clusters. To avoid this problem, FCC can be imposed using different numbers of clusters, then one can compute the fuzzy partition coefficient for each resulting clustering, which measures how well the clustering describes the data, and select the number of clusters that generates the smallest coefficient. Compared to K-means, an advantage of FCC is that it may be more flexible because it allows assigning multiple clusters to a point, thus generating more reliable results. However, FCC is more expensive in computation due to the same characteristic. There are some algorithms built on FCC, aiming at reducing the computational cost e.g., [97,98].

Fuzzy clustering has been applied to cluster catalogs e.g., [99,100], spectral data e.g., [101], and images e.g., [102]. FCC has been applied to cluster time series data e.g., [103], catalogs e.g., [104], and images e.g., [105].

3.7. Examples of Applications of Different Clustering Algorithms

In this section, we present examples of applying different clustering algorithms to the same real observational datasets in astronomy as those introduced in §2.5. Figure 10 shows the results of the clustering algorithms, including GMM, K-means, HC, DBSCAN, FCC, and self-organizing map (SOM, see §4.1). PCA is applied for visualization after the data set is clustered. The clustering algorithms are applied to generate four clusters, as suggested by the BIC calculation discussed in §3.1.

Among the six algorithms, GMM and DBSCAN yield similar results by recognizing outliers, while K-means and HC produce similar clustering. GMM, K-means, and HC identify three relatively dense, convex-shaped clusters and one scattered cluster that likely represents outliers. In contrast, DBSCAN excludes the outlier cluster from its core clusters, instead labeling the scattered points as noise, resulting in a total of five distinct clusters. Comparing GMM and DBSCAN, the overlap between corresponding clusters shows high consistency, where 83.3% of the data points are similarly clustered: the upper clusters overlap by 19.9% of the GMM upper cluster and 100.0% of the DBSCAN upper cluster; the lower right clusters by 99.3% and 96.8%, respectively; the lower left clusters by 74.4% and 99.3%, respectively (where the green and orange clusters are combined for easier comparison); the broader (i.e., outlier) clusters by 94.9% and 63.8%, respectively. Comparing K-means and HC, the overlap between corresponding clusters also shows high consistency, where 94.16% of the data points are similarly clustered: the upper clusters (pink) overlap by 91.7% of the K-means upper cluster and 85.8% of the HC upper cluster; the lower right clusters (blue) by 95.1% and 99.2%, respectively; the lower left clusters (green) by 94.6% and 93.9%, respectively; the scattered clusters on the right by 93.1% and 87.1%, respectively. The statistics may change slightly if the program is re-run, as most algorithms are not deterministic.

4. Neural Network

This section introduces two neural network algorithms: the self-organizing map and the auto-encoder. Both techniques transmit outputs between interconnected neurons that mimic the human brain, performing tasks such as clustering, dimensionality reduction, and outlier detection, or producing information that can be used as input for other algorithms.

4.1. Self-Organizing Map

Self-organizing map SOM, [106] is a neural network technique typically used for visualization by dimensionality reduction, but SOM can also be applied in clustering and outlier detection. SOM uses competitive learning, which is a form of unsupervised learning. When used for dimensionality reduction, the output is almost always two-dimensional. The number of output nodes is first defined manually, usually $k\times k$. Having more data points indicates that more output nodes, and thus a higher k, are needed. Each node is assigned a weight vector, which can be viewed as a coordinate in the input space that the node is responsible for. Therefore, the weight has the same dimensionality as the input space. Then, for each data point, SOM updates the weights of the closer nodes so the nodes become even closer to the data point. That is, the closest nodes are dragged the most, while the furthest nodes are not dragged much. The process of updating each node is iterated for the weight vectors to converge. Figure 11 shows a visualization the final result. In the end, each data point can be assigned to a winning neuron, which is the node whose weight vector is closest to the data point in the input space. Therefore, each data point can be represented by the x and y coordinates of the corresponding winning neuron. When used in clustering, the number of nodes is set to be the known number of clusters, and the data points with the same winning node belong to the same cluster. SOM performs clustering by partitioning the space, producing notably different results compared to the other clustering algorithms discussed previously, as shown in Figure 10.

Notably, SOM performs a simplification of data, not only reducing the dimensionality but also changing continuous data into discrete data. In other words, similar data points, corresponding to the same node, are presented by the same box in the map. Consequently, SOM presents the distribution of data in the 2D space clearly, but the resulting output may not be used for further data mining because most of the information is lost. Another disadvantage is that, like other neural networks, SOM is also computationally expensive.

SOM has been applied to photometric and spectroscopic data e.g., dimensionality reduction in [107,108,109,110], catalogs e.g., dimensionality reduction and clustering in [111], and light curves e.g., for clustering in [112,113].

4.2. Auto-Encoder and Variational Auto-Encoder

Auto-encoder AE, [114,115] is a neural network technique that aims to reduce a high-dimensional data set to a low-dimensional representation (a.k.a., latent space) so that the data set reconstructed from the representation is highly similar to the original data set. Therefore, AE can be used in dimensionality reduction. As shown in Figure 12, an AE network can be divided into three components: encoder, bottleneck, and decoder. The encoder consists of multiple layers, with a decreasing number of neurons in each layer, until the information reaches the bottleneck - the most compressed representation of data, with the lowest dimensionality in AE, only preserving the significant features. The bottleneck is thus the output of the encoder. If the number of neurons of the bottleneck is less than the dimensionality of the data, then the AE compresses the data [15]. The decoder works to reconstruct the data from the bottleneck. Once the encoder and decoder are trained, new data can be entered without retraining.

Similar to SOM, AE can also be used for denoising, other than dimensionality reduction. However, AE is extremely computationally expensive and may require a GPU to run. Another disadvantage of AE is that AE does not guarantee a continuous interpolation of new data because of the compactness of the latent space [15]. Therefore, variational auto-encoder VAE, [116] is introduced to solve the interpolation issue. VAE imposes a Gaussian prior, mapping each data point to a probability distribution in latent space rather than a single point. Data for reconstruction is then sampled from this distribution, mitigating issues related to interpolation and overfitting.

AE has been be applied to compress and decompress images e.g., [117,118] and catalogs e.g., [119] and to denoise time series e.g., [120]. VAE has been applied for anomaly detection of spectroscopic data e.g., [121] and time series data e.g., [122].

5. Other Applications of Unsupervised Machine Learning

This section introduces two additional applications of unsupervised machine learning: anomaly detection and symbolic regression. Anomaly detection refers to identifying outliers in a data set using various algorithmic approaches. Symbolic regression is a task where algorithms, typically genetic programming, search for one or more analytical expressions that model the data, with or without prior physical knowledge.

5.1. Anomaly Detection

Some dimensionality reduction algorithms, such as LLE and t-SNE, are good visualization tools to identify the outliers [16]. The projection often maps outliers far from the clusters, so clustering algorithms can then be applied to effectively identify the outliers e.g., [48]. In addition, some clustering algorithms can also identify the outliers. The clustering algorithm, GMM, can group the outliers into one or multiple clusters but may not classify them as outliers, as shown in Fig. Figure 10. An extension of GMM computes the log-likelihood of each sample, where a smaller score indicates a higher chance of being an outlier. The threshold for outlier detection may be set manually to achieve better results, especially when the probabilistic distribution of the data points is not Gaussian. Other clustering algorithms, such as DBSCAN and its variant HDBSCAN, can also identify the outliers, though both require certain parameters to be set manually. SOM denoising also performs outlier detection but requires a manually defined probabilistic threshold, similar to GMM.

In addition to the algorithms discussed above, Isolation Forest iForest, [123] is a dedicated anomaly detection algorithm. iForest uses random tree to randomly partition the data set. It randomly selects a feature (i.e., a dimension), chooses a random splitting threshold on that feature, and recursively partitions the data until each point is isolated in its own partition. The outliers are partitioned first because they are more isolated (i.e., sparse and located at the edges of the data), resulting in a smaller isolation path length, which is the number of partitions required to isolate a point. This suggests that, given the path length, a shorter path length suggests a higher likelihood of the point being an outlier. The path length greatly relies on the initial random partitions, so multiple trees are considered. Anomaly score is computed from the lengths obtained from multiple trees, with outliers having higher anomaly scores. The number of trees and the threshold for determining outliers are user-defined, where the threshold is often indicated by the contamination rate hyperparameter, which can significantly impact the results. Due to its structure, iForest is fast and efficient. However, iForest may not perform well when there exists interdependence between features [124]. iForest has been applied to detect the outliers in time series e.g., [124], catalogs e.g., [125], light curves e.g., [89], and latent space of neural network e.g., [126,127,128].

Figure 13 shows the exemplary results of the outlier detection using GMM, DBSCAN, SOM, and iForest. As shown, all four outlier detection algorithms identify the three major convex-shaped clusters, as discussed in §3.7, generating similar results with minor differences.

5.2. Symbolic Regression

Symbolic regression (SR) is an algorithm that computes analytical equations solely from the data, by finding the equations and their parameters simultaneously [129]. SR is based on the concept of genetic programming, which is a subfield of evolutionary computing. SR randomly initializes the first population, generating multiple configurations (i.e., equations). The symbolic expression of the equation can be visualized using tree structures with nodes and branches, each node being a symbol (e.g., ÷, log, 4.2, X, Y). SR examines the configurations and removes the less effective ones. Then, SR randomly selects and exchanges sub-configurations from two or more configurations. The evolutionary procedure iterates until SR produces a robust equation.

SR is frequently used as a supervised algorithm, where a label y is given for every data point $\overrightarrow{x}$, and SR finds the relation $y=f\left(\overrightarrow{x}\right)$. It can also be applied in an unsupervised setting, where the goal is to find a relation $f(\overrightarrow{x},y)=0$ that describes the data set [130]. There are ample applications of SR in astronomy e.g., [131,132,133,134,135], while unsupervised SR appears to be relatively unexplored.

6. Conclusion

This review discusses unsupervised machine learning algorithms used in astronomy, classifying the algorithms into three categories: dimensionality reduction, clustering, and neural networks. Anomaly detection and symbolic regression are also briefly discussed. For each algorithm, the mechanism, characteristics (e.g., advantages and disadvantages), and past applications are reviewed. Most algorithms are frequently used in astronomy, such as DBSCAN and variational auto-encoder, while some others are underutilized, such as unsupervised symbolic regression. Table 1 and Table 2 present the results of a focused search, showing the number of refereed astronomy papers in the SAO/NASA Astrophysics Data System (ADS) that applied each algorithm to different types of data from 2015 to 2025. These results offer insight into which algorithms are most widely used for specific types of analyses, reflecting current trends and preferences in the field. The rows are ordered by the sum of all algorithm applications, with higher totals ranked first.

This review also includes examples that demonstrate the results of applying these algorithms to a five-dimensional astrometry data set. Overall, unsupervised machine learning has wide application in astronomy, allowing us to analyze a large quantity of high-dimensional, unlabeled data, meeting the needs arising from technological development regarding observation.