Machine learning and Kolmogorov analysis to reveal gravitational lenses

TL;DR

This paper introduces an automated method combining Kolmogorov stochasticity and machine learning PCA to detect, classify, and catalog gravitational lenses and other astronomical objects in large datasets.

Contribution

The novel approach integrates Kolmogorov analysis with PCA for efficient detection and classification of gravitational lenses in astronomical data.

Findings

Successfully identified and classified gravitational lenses

Generated a catalog of potential lensing objects

Demonstrated high accuracy in object detection and classification

Abstract

We present an automated approach to detect and extract information from the astronomical datasets on the shapes of such objects as galaxies, star clusters and, especially, elongated ones such as the gravitational lenses. First, the Kolmogorov stochasticity parameter is used to retrieve the sub-regions that worth further attention. Then we turn to image processing and machine learning Principal Component Analysis algorithm to retrieve the sought objects and reveal the information on their morphologies. We show the capability of our automated method to identify distinct objects, including of and to classify them based on the input parameters. A catalog of possible lensing objects is retrieved as an output of the software, then their inspection is performed for the candidates that survive the filters applied.

Tables

Table 1. Table 1: List of objects: coordinates of centers, eccentricity and field numbers. 2 σ 𝜎 \sigma is used for data cut-off.

X center	Y center	Eccentricity	Field Number
306.00	289.63	0.9639	20
327.22	259.62	0.666	20
308.22	271.63	0.666	20
344.52	256.28	0.8309	21
348.22	247.62	0.666	21
338.75	252.67	0.6855	21
391.87	251.87	0.375	21
392.38	263.70	0.5264	21
458.49	315.27	0.8417	22
309.15	370.75	0.9433	28
343.40	393.57	0.7568	29
415.66	387.00	0.62	29
412.51	399.09	0.672	29
390.16	403.96	0.3564	29
458.41	340.05	0.542	30

Equations

$$F_n(x)=\{\begin{array}{ll}0,& x<X_1;\\ k/n,& X_k\le x<X_{k+1},k=1,2,\dots ,n-1;\\ 1,& X_n\le x,\end{array}$$

$$\Phi (\lambda )=\sum _{k=-\infty }^{+\infty }(-1{)}^k{e}^{-2k^2{\lambda }^2},\lambda >0,$$

$${\lambda }_{thres}={\lambda }_0+n\ast \sigma ,$$

$$\sqrt{(x_i-x_j{)}^2+(y_i-y_j{)}^2}\le \sqrt{2}.$$

$$p_i=x_i,y_i:i\in N.$$

$$\sum _{i=0}^M\sum _{j=0}^N{x}_i^m{y}_j^{n}f(x,y).$$

$$su{m}_x=\sum _{i=0}^M{x}_{i}f(x,y)$$

$$su{m}_y=\sum _{j=0}^N{y}_{i}f(x,y)$$

$$c_x=\frac{su{m}_x}{{\mu }_0},$$

$$c_y=\frac{su{m}_y}{{\mu }_0},$$

$$Co{v}_{jk}=\frac{1}{m}\sum f(x_i,y_i)(x_i^j-c_x)(y_i^k-c_y).$$

$$e=\frac{{\lambda }_1-{\lambda }_2}{{\lambda }_1}$$

$$e\ge e_{thresh}$$

$${\lambda }_2\le d_{thick}$$

Full text

††pubyear: 2019††pagerange: Machine learning and Kolmogorov analysis to reveal gravitational lenses–7

1 Introduction

The currently available various large scale sky surveys provide deep and well sampled astronomical datasets. In view of ever increasing amount of the digitized information and of the number vs their morphology of the involved sought objects, the development of efficient automatic processing of the datasets has no alternatives. Various automated methods have been described so far Alard (2006); Lenzen et al (2004); Seidel & Bartelmann (2007) to process astronomical datasets. Among them machine learning and neural network techniques Hezaveh et al (2016, 2017); Petrillo et al (2017) are becoming ever more common and powerful application means in order, first, to separate the signals of astronomical object from noise and then to classify them according to certain criteria.

The observational information on the gravitational lensing Schneider et al (1992); Straumann et al (1998); Schneider et al (2006) has become an important tool for tracing the large scale structure of the Universe, especially due to the ability for detection of the dark matter and even testing of modified gravity models, e.g. Gurzadyan & Stepanian (2018, 2019). It enables to reveal the key properties of the extragalactic objects, both of the lensed ones and of those acting as lenses, as predicted already by Zwicky Zwicky (1937). Currently morphological variety of the lens caustics, starting from twin images of the quasar SBS 0957+561 up to Einstein rings, crosses, arcs, multiple images, are discovered and interpreted, see Schneider et al (2006); Treu (2010) and references therein.

The search of the gravitational lensing evidences in the galaxy surveys includes combined use of available observational information on the lensed images and the lens itself, i.e. the spectroscopy, photometry, color, morphology. The visual inspection, however, still remains among the important steps in the recognition of elongated/distorted structures as candidates for caustics with subsequent verification by other means Wisotzki et al (2002); Frey et al (2003); Schneider et al (2006); Lopez-Caniego et al (2013); Inoue et al (2015); Mediavilla et al (2016); Nierenberg et al (2017).

Below for the first time we apply automated strategy for detection of caustics of gravitational lensing using the method of Kolmogorov stochasticity parameter (KSP) Kolmogorov (1933); Arnold (2008, 2009a, 2009b). That approach enables one to analyze signals which contain both correlated and random subsignals. Besides the study of generated sequences and modeling, that method has been already applied to the cosmic microwave background (CMB) datasets and enabled to separate e.g. the Galactic foreground, point sources (galaxies, quasars) from the cosmological signal, to analyze tiny properties of the latter Gurzadyan & Kocharyan (2008, 2009); Gurzadyan et al (2008, 2009, 2014). Among applications to datasets of quite different origin are e.g. the revealing of galaxy clusters in XMM-Newton’s X-ray data Gurzadyan et al (2011), the effect of thermal trust perturbing the trajectories of laser ranging satellites Gurzadyan et al (2013), the detection of somatic mutations in genomic sequences Gurzadyan et al (2015).

In this methodical paper we give a description of a developed three-step approach, including its application to real astronomical data. The aim of the paper is to show, first, the ability of filtering of the sub-regions that contain astronomical objects, then to identify and extract additional morphological information on those using the machine learning Principal Component Analysis algorithm.

2 The method

First, let us define the structure of the dataset under consideration. Consider an image given by 2D matrix of rectangular data, each $(x_{i},y_{i})$ pair representing a pixel with intensity $I_{ij}$.

For quantitative detection of gravitational arcs in astronomical datasets we apply 3 main steps:

Kolmogorov analysis, 2. 2.

object identification, 3. 3.

object classification.

Note that, before the first step one should gain idea on data distribution creating the histogram of intensity. That procedure is also required for Kolmogorov maps which we obtain as result of KSP analysis Kolmogorov (1933); Arnold (2008, 2009a, 2009b). This steps are needed to subtract the background from the original data, as well as to cut off the valuable data according certain criteria above given threshold. The “survived” data are analysed for object finding and then for morphological parameter obtaining.

3 Kolmogorov stochasticity parameter

Consider a real-valued sequence of numbers $\{X_{1},X_{2},\dots,X_{n}\}$ represented in increasing order. One can define two distribution functions, i.e. an empirical distribution function as given Kolmogorov (1933); Arnold (2008)

[TABLE]

and a theoretical (cumulative) distribution function as the probability $F(x)=P\{X\leq x\}.$ The difference of both distribution functions is represented by the Kolmogorov stochasticity parameter (KSP) $\lambda_{n}$ as $\lambda_{n}=\sqrt{n}\ \sup_{x}|F_{n}(x)-F(x)|.$

Kolmogorov’s theorem Kolmogorov (1933) states that \lim_{n\to\infty}P\{\lambda_{n}\leq\lambda\}=\Phi(\lambda)\, with $\Phi(0)=0$, and where

[TABLE]

so that the function $\Phi$ is independent on the theoretical distribution function $F$. The form of Kolmogorov’s function $\Phi$ determines $0.3\leq\lambda_{n}\leq 2.4$ for the KSP interval as the measure of degree of randomness of the above defined sequences Arnold (2008, 2009a, 2009b). The importance of this descriptor is that it is applicable even to sequences of few tens of length Arnold (2009a), which is not the case for most of statistical methods and is rather sensitive to the deviation from randomness. These features of the descriptor appear to be efficient at non-linear data analysis, see Atto et al (2013); Rossmanith (2013).

We then applied KSP-analysis to the observational data of a strong lensed object, namely, to the SDP.81 galaxy ALMA (2015); Tamura (2015). Our task is to find out whether KSP-method can distinguish the signal of a lensed object from that of the surrounding field, as it had enabled to separate the contribution e.g. of the Galactic disk from CMB in WMAP or Planck maps (see Gurzadyan et al (2009); Gurzadyan & Kocharyan (2009); Gurzadyan et al (2014), also for further details of the application of KSP-method).

We split the data field of size M $\times$ N into smaller sub-regions of, say m $\times$ n size, then we calculate KSPs for each sub-regions, composing so-called Kolmogorov map of the whole field. The knowledge on the original data distribution is crucial for two reasons: first, for choice of theoretical distribution function and, second, for background subtraction in the object identification step.

The detailed analysis revealed the Gaussianity of the data, hence the theoretical cumulative function was taken as Gaussian (cf.Gurzadyan & Kocharyan (2008)) and the results of the obtained values of KSP for each of the pixelized sub-regions are given in Figure 1. Before making conclusions it is worth to take a glance at the $\lambda$ distribution of 32 sub-regions (see Figure 2 ). The original data field of 672$\times$440 pixels was split into 32 sub-regions of 84$\times$110 pixels size. Therefore, in the Figure we show the histogram made of 32 $\lambda$ values, having mean value of around 1.9. One can clearly see that the majority of $\lambda$ values are between 0.5 and 2.2, which confirms the correctness of our assumption regarding the data Gaussianity. Higher values are due to the more regular structures in the sub-regions. The sub-regions that contain parts of lensing arcs have anomalously high $\lambda$.

The further steps strongly depend on the data type: if the original image is sparse, it is appropriate to take the modal value (although, in this case the mean is close to modal) for $\lambda$ as background, otherwise, if the image is full of objects, the mean value for $\lambda$ might be considered. We filter KSP map with cut-off value defined as

[TABLE]

where $\lambda_{0}$ is either mean or modal value of the entire field, $\sigma$ is the standard deviation, and finally n multiplier indicates how many $\sigma$-s we want to cut above the mean or modal value. This multiplier should rather be decided empirically.

Due to the data spread we get large value of standard deviation, therefore $n=1$ is a good choice for multiplier. This results in 6 sub-regions to pass our filter. By next step we are going to identify objects in those regions.

4 Object identification

This section we devote to description of the object finding algorithm. Let us first define the object. It is a set of one-connected pixels that are isolated from other sets of connected pixels. And as far as we are looking for astronomical object composing pixels, those pixels should have somewhat higher values of intensity. Therefore using our knowledge about original data (distribution, mean or modal value and standard deviation) we apply cut-off technique and maintain only pixels that survive the filter (survived data). Having a list of those pixels p($x_{i},y_{i}$) we define Moore neighborhood N(p($x_{i},y_{i}$)) for each of them, that is the set of all pixels that are orthogonally or diagonally-adjacent to the region of interest and the region of interest itself may or may not be considered part of the Moore neighborhood Moore (1964). In other words, for the given pixel with ($x_{i},y_{i}$) coordinates, the neighbors are those with coordinates obeying the following inequality

[TABLE]

There are different approaches to check the path connectivity between two pixels Soille (2003); Lenzen et al (2004). However we propose the simple algorithm that we came up with.

First of all we construct the lists of neighbors for each pixel. As far as we search those pixels from the survived data, the boundary pixels of objects do not contain the full set of neighbors. After this step, we take an auxiliary list A that, in the end is supposed to contain the sets of pixels for different objects (it is set of those sets), and compare its content with every N(p).

$a)\,$ If $A\cap N(p)=\{p_{k}...p_{m}\}$ , this means that current pixel belongs to object $\to$ add $A\setminus N(p)$ into A;

$b)\,$ Otherwise, if $A\cap N(p)=\emptyset$, $\to$ iterate over elements of A, if finds intersecting lists, $\to$ go to $a)\,$, completes, then breaks, otherwise if it iterates till the end, just adds a new set into A with the current N(p) elements and then breaks.

This procedure results with the set of isolated objects for each sub-regions.

5 Object classification

Object classification is the final and the most important step of this algorithm. Up to now we described our approach to the automated arc searching using the KSP-method that helps saving computational efforts by filtering the fields with valuable information. Then, we developed an object finding algorithm and now we are going to give detailed description of the object classification algorithm.

First of all, let us recall that, the gravitational lensing arcs are elongated objects - particularly arcs of a circle originated due to the lens distribution along the line-of-sight. Hence, our task is to search for that kind of objects, sometimes as tiny as PSF. Some galactic spiral arms are similar to arcs, so can be mutually confused, if not the typically significant thickness of arms.

Let the isolated object be composed of set of pixels

[TABLE]

The mathematical moment of m+n order is defined as

[TABLE]

In order to calculate the object’s total intensity and the center we calculate 0-th and 1-th order moments, correspondingly,

[TABLE]

and

[TABLE]

and then, we get the center of the object as follows

[TABLE]

[TABLE]

where $\mu_{0}$ is the total intensity of the object. The second order moment is defined as $2\times 2$ matrix with components

[TABLE]

This matrix is called covariance matrix. It is a symmetric matrix having variance values for x and y independent variables in its diagonal. While the variance is defined for a set of variables and refers to the spread of data around its mean, the covariance refers to the measure of the directional relationship between two random variables. Hence it can be used to investigate the properties of isolated objects, whether elongated or more regular. To this end we calculate the eigenvalues of matrix $Cov_{jk}$ and eigenvectors. Whereas the eigenvectors define the principal components of the data, the eigenvalues are the scales along those principal components.

A $2\times 2$ matrix can have at most two eigenvalues, say $\lambda_{1}$ and $\lambda_{2}$, and if $\lambda_{1}$ $>$ $\lambda_{2}$, then

[TABLE]

is the eccentricity of the object. Our software includes threshold parameter for eccentricity and thickness of the elongated object. So, depending on data we study, those threshold values vary and might be induced empirically. Hence, at the end, our software outputs a list of objects with a diagnosis about object type. Of course the results depend on our choice of $e_{thresh}$ and $d_{thick}$. The objects for which

[TABLE]

and

[TABLE]

only appear in the final list.

6 Results

Below we introduce a table that contains a list of isolated objects. We set the threshold value for eccentricity 0.35, therefore in the table we have objects with small eccentricities. The reason we keep objects with small eccentricities is that, in the first step we split the whole field into smaller sub-regions, and as it has been shown in Table 1, the lensing structure is separated between parts. This assumes an extra work of grouping the parts of the same structure together, but at the same time we make profit when with Kolmogorov analysis we avoid performing the search of structures in the entire field. Indeed, we just look for objects in the fields that survive $\lambda$ filtering.

From the Table 1 one can see that the first 3 objects belong to the same field, and from the first glance those are in the same big arc. However attentive inspection of Fig. 3, in particular of sub-region 20, makes clear that along the big arc there is variation of pixel intensities, so that some part does not pass the filter and the visibly continuous arc splits to several disjoint arcs of smaller size. In Fig. 4 we represent the results of our software for different cut-off values of intensities, particularly, Fig. 4(a) corresponds to 1-$\sigma$, Fig. 4(b) to 1.5-$\sigma$, Fig. 4(c) to 2-$\sigma$ and Fig. 4(d) to 3-$\sigma$.

Note that, the results of this method depend on the quality of image. Due to the sensitivity of KSP-method, even the hidden from eye signals can be revealed. Even if the background of field is quasi-uniform at high level, small perturbations can be detected by KSP-method. An efficient Object identification and Object classification are matter of careful understanding of the data available. Indeed, for those steps it is crucial to know the data distribution, as well as its mean or modal values, and to perform a proper cut-off to avoid loosing valuable information.

Additionally, the ability of the suggested method to reveal "low visibility" structures can be seen from the following example. We simulated images of different statistical significance. Namely, even if the object pixels’ intensities are only 0.5-$\sigma$ above the background level, this method is still able to identify elongated objects (Fig. 5). Indeed, it is hardly possible to notice structures in the initial simulated image (left plot), whereas the right image presents the results retrieved by the algorithm described above. Of course, as one might expect, other structures may be retrieved as well as shown in the right plot and their possible association to the lensing structure candidates should be investigated additionally.

7 Conclusions

We advanced an automated method of search for isolated objects in astronomical datasets, that is, extraction of valuable information based on statistical properties of data and Moore’s neighborhood algorithm.

First, the split of pixelized regions containing the gravitational lens images showed statistically notable difference regarding the Kolmogorov function with respect to their average value in the surrounding sky. Then, we extracted morphological information on the extracted objects applying Principal Component Analysis strategy. This enabled us to classify objects among elongated or regular structures. The efficiency of the method is illustrated via simulation of low significance structures.

Observational surveys offer huge amount of data and development of any new tool for automated search for given category of objects can be important. Particularly, the Kolmogorov analysis in addition to the problems mentioned in the Introduction can be applied to test the isotropy of sky distribution of gamma-ray bursts. That problem has been studied by various methods, not always with identical conclusions (see Ruggeri & Capozziello (2016); Andrade et al (2019); Ripa & Shafieloo (2019) and references therein), and is of remarkable cosmological importance. Image processing and machine learning algorithms are currently becoming conventional tools for astronomical datasets, to extract ever more refined information. In this paper we considered an example of combination of statistical methods and machine learning algorithms to reveal certain structures in astronomical datasets.

8 Acknowledgments

We acknowledge the use of data from http://almascience.nrao.edu/aq/.