LaTeX Deep Learning application for stellar parameters determination:
I- Constraining the hyperparameters

Machine learning (ML) applications have been used extensively in astronomy over the last decade ([6]). This is mainly due to the large amount of data that are recovered from space and ground-based observatories. There is therefore a need to analyze this data in an automated way. Statistical approaches, dimensionality reduction, wavelet decomposition, ML, and Deep Learning (DL) are all examples of the attempts that were performed in order to derive more accurate stellar parameters such as the effective temperature ($T_{\mathrm{eff}}$), surface gravity ($\log g$), projected equatorial rotational velocity ($v_{e}\sin i$), and metallicity ($[M/H]$) using stellar spectra in different wavelength ranges ([26, 46, 47, 59, 62, 4, 32, 15, 21, 38, 19, 42, 43]). DL is a machine learning method based on deep Artificial Neural Networks (ANN) that does not usually require a specific statistical algorithm to predict a solution but it is rather learned by experience and thus require a very large dataset ([65]) for training in order to perform properly.

An overview of the automated techniques used in stellar parameters determination can be found in [32]. We will mention some of the recent studies that involve ML/DL. The increase of the computational power and the large availability of predefined optimized ML packages (in e.g. Python, C++, and R) have allowed astronomers to shift from classical techniques to ML when using large data. One of the first trials to derive the stellar parameters using neural networks was done by [5]. This work demonstrated that networks can give accurate spectral type classifications across the spectral types range B2-M7.

[14] presented an ANN architecture that learns the function which can relate the stellar parameters to the input spectra. They obtained residuals in the derivation of the metallicity below 0.1 dex for stars with Gaia magnitude $\mathrm{G}_{rvs}<12$ mag, which accounts for a number in the order of four million stars to be observed by the Radial Velocity Spectrograph of the Gaia satellite¹¹1The limited magnitude of the RVS is around 15.5 mag ([12]).. [48] used a ML algorithm to measure the mean longitudinal magnetic field in stars from polarized spectra of high resolution. They found a considerable improvement of the results, allowing to estimate the errors associated to the measurements of stellar magnetic fields at different noise levels.

[44] developed, and applied a Convolutional Neural Network (CNN) architecture using multi-task learning to search for and characterize strong HI Ly$\alpha$ absorption in quasar spectra. [15] applied a deep neural network architecture to analyse both SDSS-III APOGEE DR13 and synthetic stellar spectra. This work demonstrated that the stellar parameters are determined with similar precision and accuracy as the APOGEE pipeline.

[55] introduced an automated approach for the classification of stellar spectra in the optical region using CNN. They also showed that deep learning methods with a larger number of layers allow the use of finer details in the spectrum which results in improved accuracy and better generalisation with respect to traditional ML techniques.

[59] introduced a deep-learning method, SPCANet, that derived $T_{\mathrm{eff}}$ and $\log g$ and 13 chemical abundances for LAMOST Medium-Resolution Survey (MRS) data. These authors found abundance precision up to 0.19 dex for spectra with Signal-to-Noise ratio (SNR) down to $\sim$10. The results of SPCANet are consistent with those from other surveys such as APOGEE, GALAH and RAVE, and are also validated with the previous literature values including clusters and field stars. [26] derived the atmospheric parameters and abundances of different species for 420,165 RAVE spectra. They showed that CNN-based methods provide a powerful way to combine spectroscopic, photometric, and astrometric data without the need to apply any priors in the form of stellar evolutionary models.

More recently, [35] introduced a multi-layer CNN to forecast solar flare events probability occurrence of M and X classes. [10] introduced an AGN recognition method based on Deep Neural Network. [1] used machine learning ML methods to generate model SEDs and fit sparse observations of low-luminosity active galactic nuclei. [50, 49] used CNNs and different ANN networks to estimate emission-line parameters and line ratios present in different filters of SITELLE spectrometer. [13] used DL combined with k-Nearest Neighbour and Decision Tree Regression algorithms to compare the accuracy of the predicted photometric redshifts of newly detected sources. [41] applied the ThetaRay Artificial Intelligence algorithms to 10 803 light curves of threshold crossing events and uncovered 39 new exoplanetary candidates targets. [8] reached a classification accuracy of 88 percent while investigating the use of a CNN for automated merger classification. [16] used an assisted inversion techniques based on CNN for solar Stokes profile inversions. In the context of Classification of galactic morphologies, [17] used a ML generative adversarial networks to convert ground-based Subaru Telescope blurred images into quasi Hubble Space Telescope images. [18] presented StelNet, a Deep Neural Network trained on stellar evolutionary tracks that quickly and accurately predicts mass and age from absolute luminosity and effective temperature for stars of solar metallicity.

In this manuscript, we present both a new method to derive stellar atmospheric parameters, and we also demonstrate the effect of each of the CNN parameters (such as the choice of the optimizers, loss function, activation function, $\ldots$) on the accuracy of the results. We will provide the procedure that can be followed in order to find the most appropriate configuration independently of the architecture of the CNN. This is intended as the first in a series of papers that will help the astronomical community to understand the effect on the accuracy of the prediction from most of the parameters and the architecture of the network. CNN parameters are numerous and to find the optimal ones is a very hard task. To do so, we trained the CNNs with different configurations of the parameters using purely synthetic spectra for the 3 steps of training, cross-validation (hereafter called validation), and testing. Using synthetic spectra, we have access to the true parameters during our tests. Noisy spectra are tested in order to mimic observations.

We have limited our work to a specific type of objects, A stars, because as mentioned previously the purpose is not to show how well we can derive the labeled stellar parameters but what is the effect of specific parameters on stellar spectra analysis. By applying our models to A stars, we use previous results ([19, 32]) as a reference for the expected accuracy of the derived stellar parameters. In the same way, the wavelength range and the resolving power are chosen to be representative of values used by most available instruments. Once the calibration of the hyperparameters was performed, we have tested our optimal network configurations on a set of FGK stars in Sec. 6, using the wavelength range of [42].

The training, validation, and test data are explained in Sec. 2. Section  3 discusses the data preparation previous to training. The neural network construction and the parameters selection is explained in Sec. 4. Results are summarized in Sec. 5. The application of the optimal networks to FGK stars is performed in Sec. 6. Discussion and conclusion are gathered in Sec. 7.

Our learning, or training databases (TDB) are constructed from synthetic spectra for stars having effective temperature between 7 000 and 10 000 K, and the wavelength range of 4450 Å to 5000 Å. This range was selected because it is in the visible domain and contains metallic and Balmer lines sensitive to all stellar parameters ($T_{\mathrm{eff}}$, $\log g$, $[M/H]$, $v_{e}\sin i$), especially for the spectra types selected in this work. This region is also insensitive to microturbulent velocity which was adopted to be $\xi_{t}$=2 km/s based on the work of [19, 20]. Surface gravity, $\log g$, is selected to be in the range of 2.0–5.0 dex. Projected rotational velocity, $v_{e}\sin i$, is calculated between 0–300 $\mbox{km}\,\mbox{s}^{-1}$. The Metallicity, $[M/H]$, is in the range of -1.5 and +1.5 dex. Table  1 displays the range of all stellar parameters. These spectra are used for both the training and the validation phases. Approximately 55 000 noise free synthetic spectra were calculated using a random selection of the stellar parameters in the range of Tab. 1. These spectra are used instead of the observations (Test data without noise). Gaussian SNR, ranging between 5 and 300, were added to these test spectra in order to check the accuracy of the technique on noisy data (Test data with noise).

Details for the calculations of the synthetic spectra can be found in [19] or [32]. In summary, 1D plane–parallel model atmospheres were calculated using ATLAS9 ([34]). These models are in local thermodynamic equilibrium (LTE) and in hydrostatic and radiative equilibrium. We have used the new opacity distribution function in the calculations ([9]) as well as a mixing length parameter of 0.5 for 7000 K$\leq$$T_{\mathrm{eff}}$$\leq$8500 K, and 1.25 for $T_{\mathrm{eff}}$$\leq$7000 K ([57]).

We have used [30] SYNSPEC48 synthetic spectra code to calculate all normalized spectra. The adopted line lists are detailed in [19]. This list is mainly compiled using the data from Kurucz gfhyperall.dat²²2http://kurucz.harvard.edu, VALD³³3http://www.astro.uu.se/$\sim$vald/php/vald.php, and the NIST⁴⁴4http://physics.nist.gov databases.

Finally, the resolving power is simulated to R=60 000. This value falls in the range between low and high resolution spectrographs. The technique that will be shown in the next sections can be used for any resolution. The construction and the size of the TDB will be discussed in Sec. 5. The use of synthetic spectra in ML to constrain the stellar parameters has shown to suffer from the so-called synthetic gap ([15, 46]). This gap refers to the differences in feature distributions between synthetic and observed data. We have decided to limit our work to synthetic data for 2 reasons, first we would like to remove the hassle of the data preparation steps (data reduction, flux calibration, flux normalization, radial velocity correction $\ldots$), and second because our intention is to find the strategy and technique that should be adopted in ML for deriving stellar parameters.

We are working on a future paper that deals with the architecture of the network as well as the choice of the kernel sizes and the number of neurons. Combining the best strategy to constrain the hyperparameters (this manuscript) as well as the most optimal architecture (future studies) will allow us to use a combination of synthetic and observational data in our training database. Having well known stellar parameters, these observational data will allow us to remove/minimize the synthetic gap and better constrain the stellar parameters.

The TDB contains $N_{\mathrm{spectra}}$ spectra that span the wavelength range of 4450–5000 Å. Having a wavelength step of 0.05 Å, this results in $N_{\lambda}=$10800 flux points per spectrum. The TDB can then be represented by a matrix M of size $N_{\mathrm{spectra}}\times N_{\lambda}$. A color map of a subsample of M is displayed in Fig.1. Although the synthetic spectra are normalized, some wavelength points could have fluxes larger than unity. This is due to the noise that is incorporated during the so-called data augmentation procedure which will be explained in Sec. 4.1.1.

Training a CNN using the M matrix is time consuming, especially if one should use a larger wavelength range or a higher resolution. For that reason we have applied a dimensionality reduction technique i.e., Principal Component Analysis (PCA hereafter), in order to reduce the size of the training TDB as well as the size of the validation, test, and noisy synthetic data. Although this step is optional, we recommend its use whenever the data can be represented by a small number of coefficients. The PCA can reduce the size of each spectrum from $N_{\lambda}$ to $n_{k}$. The choice of $n_{k}$ depends on the many parameters, the size of the database, the wavelength range and the shape of the spectra lines. As a first step, we need to find the Principal Components, and to do so, we proceed as follows:

The matrix M is averaged along the $N_{\mathrm{spectra}}$-axis and the result is stored in a vector $\bar{M}$. Then, we calculate the eigenvectors $\textbf{e}_{k}(\lambda)$ of the variance-covariance matrix C defined as

where the superscript ”T” stands for the transpose operator. C has a dimension of $N_{\lambda}\times N_{\lambda}$. Sorting the eigenvectors of the variance-covariance matrix in decreasing magnitude will results in the ”Principal Components”. Each spectrum of M is then projected on these Principal Components in order to find its corresponding coefficient $p_{jk}$ defined as

The choice of the number of coefficient is regulated by the reconstructed error as detailed in [42]:

We have opted to a value for $n_{k}$ that reduces the mean reconstructed error to a value ¡0.5%. As an example, using a database of 25 000 spectra with stellar parameters ranging randomly between the values in Tab. 1 requires less than 7 coefficients to reach an accuracy ¡1%, and a value of $n_{k}$=17 to reach a 0.5% error as shown in Fig. 2. This technique has shown its efficiency when applied to synthetic and/or real observational data with $T_{\mathrm{eff}}$$>$ 4000 K (see [19, 42, 43], for more details).

Applying the same procedure to all our TDB and taking the maximum value to be used for all, we have adopted a constant value for $n_{k}$=50. This value takes into account all the databases that will be dealt with in this work, especially that some will be data augmented as will be explained in Sec. 4.1.1. This means that instead of training a matrix having a dimension of $N_{\mathrm{spectra}}\times N_{\lambda}$, we are using one with dimension of $N_{\mathrm{spectra}}\times n_{k}$, with $n_{k}\ll N_{\lambda}$. In that case, our new data consist of a matrix containing the coefficients that are calculated by projecting the spectra on the $n_{k}$ eigenvectors.

This projection procedure over the Principal Components is then applied to the validation, test, and noisy spectra datasets.

This section begins with a brief description of supervised⁵⁵5Supervised learning refers to algorithms that calculate a predictive model using data points with known labels/outcomes. learning. Given a data set $(X,Y)$, the goal is to find a function $f$ such that $f(X)$ is as ”close” as possible to $Y$. For example, $Y$ could be the effective temperatures or surface gravity and $X$ the corresponding spectra. This ”closeness” is typically measured by defining a loss function $L(f(X),Y)$ that measures the difference between the predicted and actual values. Therefore the goal of the learning process is to find $f$ that minimizes $L$ for a given dataset $(X,Y)$. Ultimately, the success of any learning method is assessed by how well it generalizes. In other words, once the optimal $f$ is found for the training set $(X,Y)$, and given another data set $(U,V)$, how close is $f(U)$ to $V$.

One of the most successful methods in tackling this kind of problems is Artificial Neural Networks (ANN), a subset of ML. As the name suggests, an ANN is a set of connected building blocks called neurons which are meant to mimic the operations of biological neurons ([2, 60]). Different kinds of ANNs can be built by varying the number of, connections between, and operations of individual neurons. The operations performed by these neurons depend on a number of parameters, called weights and some nonlinear function, called the activation. At a high level, an ANN is just the function $f$ that was described earlier. Since the network architecture is chosen at the start, finding the optimal $f$ boils down to finding the optimal weight parameters that minimize the cost function $L$.

Regardless of the type of ANN used, the process of finding the optimal weights is more or less the same, and works as follows. After the network architecture is chosen, the weights are initialized, then a variant of gradient descent is applied to the training data. Gradient descent changes the parameters iteratively, at a certain rate proportional to their gradient, until the loss value is sufficiently small ([52]). The proportionality constant is called the learning rate. While this process is well known, there is to date no clear prescription for choice of the different components . The main difficulty arises from the fact that the loss function contains multiple minima with different generalization properties. In other words, not all minima of the loss function are equal in terms of generalization. Which minimum is reached at the end of training phase depends on the initial values chosen for the weights, the optimization algorithm used, including the learning rate and the training dataset ([63]). In the absence of clear theoretical prescriptions for the components one has to rely on experience and best practices ([7]).

One popular type of ANN is the feedforward network, where neurons are organized in layers, with the outputs of each layer fully connected to the inputs of the next. By increasing the number of layers (whence the ”deep” in ”deep learning”) many types of data can be modelled to a high degree of accuracy. Fully connected ANN, however, have some shortcomings, such as the large number of parameters, slow convergence, overfitting, and most importantly, failure to detect local patterns. Almost all the aforementioned shortcomings are solved by using convolution layers.

The effect of each CNN parameter on the accuracy of the stellar parameters has been tested. To do so, we have used the same CNN with the same parameters for all our tests while changing only the concerned one at each time. For example, to find the best epoch numbers, we fix the Activation function, the optimizer, the number of Batches, the Dropout percentage, the loss function and the kernel initialiser while iterating on the number of epochs. The same parameters are used again for finding the optimal dropout percentage and so on. The fixed values used in these calculations are the he_normal for the kernel initialiser, the mean squared error for the loss function, the ”ADAM” optimizer, the relu activation function, 50% of Dropout, 64 Batches. These tests are performed with epochs of 100, 500, 1000, 2000, 3000, 4000 and 5000. In all tests, the distribution of Training and Validation are 80% and 20% respectively.

The results will be a combination of Test errors spanning over different number of epochs for each stellar parameter and CNN configuration. The variation with the number of epochs ensures that the trends are real and not due to local minima as a result of the low number of iterations. The Tests are a collection of 110 000 synthetic spectra, half of them without noise and half with random noise as introduced in Sec. 2.

To better visualise the results and to have a better conclusion about the optimal configurations, we display in Figs. 5 to 8 the relative error of the observations. These errors are calculated by dividing the values by the maximum observation standard deviation in all configurations (ie. including all epochs simulations). This will allow us to target the minimum values and pinpoint the best parameters.

In what follows, we show the results that were performed using a training dataset of  40 000 randomly generated synthetic spectra in the ranges of Tab. 1. In Sec. 5.5, we discuss the effect of using a small or a large training database and the effect of using or not Data Augmentation.

In order to verify how universal the results are, we checked that the optimization of the code is not dependent on wavelength and/or spectral type, we also tested the procedure on FGK stars. To do that, we have calculated a TDB specific for FGK stars using the parameters displayed in Tab. 4. The wavelength range was selected to coincide with the one of [42]. This range is sensitive to all the concerned stellar parameters.

A database of 50 000 random synthetic spectra with known stellar labels is used in the training. About 20 000 Test data, with and without noise, were calculated in the same range of Tab. 4 to be used for verification. The optimal NNs that were introduced in Sec. 5 were used again, as a proof of concept, for the FGK TDB. The results are displayed in Tab. 5 for Training, Validation, and tests.

Because of the low rotational velocities for FGK stars ($v_{e}\sin i$$<$ 100 $\mbox{km}\,\mbox{s}^{-1}$), the results are more accurate. That is not surprising because $v_{e}\sin i$ drastically affects the shape of the lines as in A stars. The derived errors on the stellar parameters are found to be 82 K, 0.07 dex, 0.90 $\mbox{km}\,\mbox{s}^{-1}$and 0.06 dex for $T_{\mathrm{eff}}$, $\log g$, $v_{e}\sin i$and $[M/H]$, respectively (Tab. 5). These results are very promising, but we should be aware of the complications that would arise when using real observations, especially in the case of the cool M stars. These stars have been analyzed in the context of exoplanets search ([54, 46]) and show complications in their spectra mainly related to the continuum normalization. Adapting the data preparation and the CNN will be inevitable in order to take into account these effects. These results also show that when deriving the stellar parameters for specific spectral types, the wavelength region should be selected according to these spectral lines/bands most sensitive to the variations of the parameters one seeks.

The purpose of this work is not only to find the best tool for the accurate prediction of parameters but also to show the steps that should be taken in order to reach the optimal selection of the CNN parameters. Often scientists use DL as a black box without explaining the choice of the parameters and/or architecture. In this manuscript, we have explained the reason for selecting specific hyperparameters while emphasizing the pedagogical approach. To have a more effective tool, one should change the architecture of the model. The architecture of the model depends on the type and range of the input. In this work we have fixed the architecture and iterated on the hyperparameters only.

Sections 5.1 to 5.4 show that for each stellar parameter, the setup of the network should be changed. This means that for a specific network and a specific stellar parameter, a study should be made to find the optimal configuration of hyperparameters. This is due to the contribution of the specific stellar parameter on the shape of the input spectrum. Using the PCA decomposition, we have reduced the size of the input parameters to only 50 points per spectrum while keeping more than 99.5% of the information. This is recommended in case of large databases and wide wavelength range and could avoid the use of extra pooling layers in the network. This projection technique is not only applicable for AFGK stars but can also be used for cooler stars. Further, ([29, 43]), [53] have applied a projection pursuit regression model based on the independent component analysis compression coefficients to derive $T_{\mathrm{eff}}$, $\log g$, and $[M/H]$ of M-type stars.

Although the CNN architecture was not optimized, we were able, using a strategy of finding the best hyperparameters, to reach a level of accuracy that is comparable to other adopted techniques. In fact, we found for A stars, an average accuracy of 0.08 dex for $\log g$, 0.07 dex for $[M/H]$, 3.90 $\mbox{km}\,\mbox{s}^{-1}$ for $v_{e}\sin i$, and 127 K for $T_{\mathrm{eff}}$. In the case of stars with $v_{e}\sin i$ less than 100 $\mbox{km}\,\mbox{s}^{-1}$, we found the accuracy to be 90 K, 0.06 dex, 0.06 dex and 2.0 $\mbox{km}\,\mbox{s}^{-1}$, for $T_{\mathrm{eff}}$, $\log g$, $[M/H]$ and $v_{e}\sin i$, respectively. These accuracy values are signal to noise dependant and reduce as long as the signal to noise increases. Extrapolating the technique to FGK stars also shows that the same network could be applied to different spectral types and different wavelength ranges.

The technique that we followed in this paper could be transferable to any classification problem that involves neural network. In the future we plan to develop a strategy to find the best CNN architecture depending on the input data and the type of the predicted parameters. Once the architecture and the configuration of the parameters is settled, we will be testing the procedure on observational spectra as we did in [42], [43], [19] and [32]. Using only observational data or a combination of synthetic spectra and real observations with well known parameters will allow us to constrain the derived stellar labels while minimizing the critical synthetic gap ([15]). One more criterion that should be taken into account is when applying this technique to real observations, thorough data preparation work should be done to take into account the characteristics of each spectral type (e.g., continuum normalization in M and giant stars, and low number of lines in hot stars.).