On the use of machine learning algorithms in the measurement of stellar magnetic fields
J. C. Ramirez-Velez C. Yañez-Marquez J. P. Cordova-Barbosa




Instituto de Astronomía - Universidad Nacional Autónoma de México, 
           Apdo. Postal 877, 22860, Ensenada B.C.,  México,  jramirez@astro.unam.mx
 Laboratorio de Cómputo Inteligente - Instituto Politécnico Nacional,  Av. Juan de Dios Bátiz s/n, CP 07738, CDMX,
México. 
 CUCEA, Universidad de Guadalajara, Periférico Norte 799, Los Belenes,  Zapopan Jalisco, 45100, México.






jramirez@astro.unam.mx



Regression methods based in Machine Learning Algorithms (MLA) have become an important tool for
data analysis in many different disciplines.

In this work, we use MLA in an astrophysical context; our goal is to measure the mean longitudinal magnetic field in stars
(H_ eff) from polarized spectra of high resolution, through the inversion of the so-called multi-line profiles.

Using synthetic data, we tested the performance of our technique considering different noise levels:
In an ideal scenario of noise-free multi-line profiles, the inversion results are excellent; however, the
accuracy of the inversions  diminish considerably when noise is taken into account. In consequence, we propose a
data pre-process in order to reduce the noise impact, which consists in a  denoising profile process combined  
with an iterative inversion methodology.

Applying this data pre-process, we have found a considerable improvement of the inversions results, 
allowing to estimate the errors associated to the measurements of stellar  magnetic fields at different noise levels.

We have successfully  applied our data analysis technique to two different stars, attaining by first time the measurement of
H_ eff from multi-line profiles beyond the condition of line autosimilarity  assumed by other techniques. 




Machine learning algorithms applied to the measurement of H_ eff.
Ramírez Vélez  et al.  
On the use of machine learning algorithms in the measurement of stellar magnetic fields
     
     J.C. Ramírez Vélez1
 C. Yáñez Márquez2
 J.P. Córdova Barbosa3 



    
=======================================================================================





§ INTRODUCTION


Magnetic fields are a key ingredient in the study of the stellar evolution 
and yet many of their aspects 
are still unkown, among  other reasons, due to the difficulty to 
detect, measure and map (whenever this is possible) 
surface magnetic fields in  stars  other than the Sun. 

Nonetheless, a lot of progress has been made in the recent years in the study
and observation of stellar magnetic fields, and  some facts are now 
well-stablished. For example, non chemically peculiar massive stars
can host stable strong dipolar configurations and the mean 
longitudinal magnetic  fields (H_ eff) for 
these stars vary from tens to thousands of gauss
<cit.>; 
or that the solar type stars during the main sequence phase, 
harbor more complex magnetic geometries and exhibit
much weaker effective fields (H_ eff), by at least one order of magnitude, 
namely, in the order of tens of gauss or less <cit.>.
These two examples are of interest for this work from a data analysis point of view, 
because their results were obtained using the so-called multi-lines 
approach. This means that to retrieve a signal in the circular polarised 
Stokes parameter (V),  in order to detect and to subsequently measure the magnetic field, 
it is necessary to co-add as many as possible individual spectral lines. The line 
addition permits to boost the signal-to-ratio value in the  resulting “mean” 
circular profile. Otherwise, the individual V profiles are in general 
systematically bellow the noise level.


The underlying idea of grouping many lines toghether into a single one is that the
noise addition is incoherent while the  addition of the V profiles is coherent 
<cit.>, i.e., the more individual profiles are co-added,  the better 
signal to noise ratio is retrieved in the mean profile. 
Typically, an ad-hoc selection of unblended lines is performed establishing a linear 
relation between then mean circular profile and the mean longitudinal magnetic field.


The  number of lines to combine depends on the spectral type of the star and on the 
instrumental resolution. In the case of modern high resolution spectroplarimeters such 
as HARPS, ESPaDOnS or NARVAL, the total number can reach up to 
ten thousand lines, or more,   for late-type stars.

Different techniques have been proposed to perform the line addition, and they can be
roughly divided into two groups: 1) the line addition is based on the 
assumption that all the  lines are autosimilars <cit.> 
or apply deblending process to extend the range of applicability of the technique
<cit.> and 2)  the 
line autosimilarity assumption is not necessar <cit.>.
The mean profiles obtained in the former approach are called  Least-Square-Deconvolved 
(LSD) profiles after the work of <cit.>,   and they require the use of unblended 
lines, while the mean profiles issued from the latter approach are called Multi-Zeeman-Signatures (MZS) 
after  the work of <cit.>, and in this case the mean profiles are established
including blended and unblended lines.


For the LSD profiles, the measurement of the strength of the mean longitudinal
magnetic fields can be done through the weak field approximation 
<cit.>, or alternatively, with the centre-of-gravity 
method <cit.> through the first-order 
moment of the mean circular profile <cit.>. 
The employment of the WFA and of the COG methods  are of great 
utility because they allow to measure the effective magnetic field withouth having 
to recourse to radiative transfer calculations. Nowadays, the COG method is the technique
used by excellence to determine H_ eff from the LSD profiles.
Nevertheless, both approaches have limitations due
to assumption of self-similarity of the lines, 
having to fulfill simoultaneously both contidions: 
the use of unblended lines and the measurement of
weak magnetic fields. The range of application in the data analysis
is thus difficult to determine precisely. 
For example, concerning the second condition related to the strength of the magnetic fields, 
<cit.> stated a maximum of 2 kG as limit for the analysis of circular
polarised LSD profiles; however, as noticed by
<cit.> the limit of the WFA of local profiles, i.e. to reach
the Zeeman saturation regime in individual lines,  is around one kG
for lines in the optical wavelength range. In any case, to bypass the inherent limitations  
associated to the use of the LSD profiles and to disposse of a more general approach, 
in this work we continue developping a method based in the analysis of the MZS profiles.


In a previous work, <cit.>, we presented a
new inversion code for the measurement of stellar magnetic fields using a 
complete radiative transfer approach. We used as case of study  a solar type 
star (T_ eff = 5 500 K), to show that the MZSs can be used to 
correctly retrieve the stellar longitudinal magnetic field in cool stars. 
The proposed methodology illustrates that it is feasible to infer H_ eff, 
in the context of the multi-lines approach, beyond the assumption of 
lines autosimiliraity, i.e, beyond the employment of the COG method (or the WFA).

The inversion code used in  was based on a look-up table, i.e, 
the higher the number of stellar MZSs contained in the table the better were the results
obtained. This fact makes of our technique a very CPU-time consuming method. 
In this work we incorporate the use of machine learning algorithms to remarkably 
decrease the number of synthetic stellar spectra required to properly determine 
H_ eff.

In section II, we show the performance of different machine learning algorithms for the 
inversion of ideal noiseless MZSs, to subsequently show that the effectivity of the
regression alrogrithms is strongly affected when noise is included, even for very small 
noise levels. In section III we incorporate pre-process data analysis  to reduce 
the noise effect when infering H_ eff. In section IV we apply our
method to two real cases and compare our results to previous ones published in others 
studies, to finally drew some general conclusions in section V.



§ THE MACHINE LEARNING STRATEGY


In the stellar pyhisics domain, the number of published studies using 
analysis methods based in machine learning algorithms has been growing,
specially in the last years, either for classification of stars
<cit.>
or in the search of exoplanets <cit.>, or in the 
determination of stellar fundamental parameters 
<cit.>, among many others.
In this work we will use the machine learning approach in the 
context of the study of stellar magnetic fields.



As previously mentioned, in  we performed inversions of MZSs
using a code based in a look-up table strategy, meaning that given a MZSs to invert 
we find as solution the closest-norm MZSs in the table.  The MZSs 
were constructed 
adopting the eccentric tilted dipole model <cit.>. In this model, we 
considered as free: 1) all the parameters that describe the magnetic geometry of 
the dipolar configuration (the three Eulerian angles α, β, γ, and the 
inclination angle (i) of the stellar rotation axis respect to the line of sigth),  
2) the position of the magnetic dipole inside of the star given by two coordinates (X_2,X_3),
and 3) the strength of the dipolar moment (m), 
leading to a total of 7 free parameters. It turned out that in order to obtain good
enough inversion results, it was required to calcute 7 500 stellar spectra 
models. For each synthetic spectrum we randomly varied the 7 free parameters,
covering a spectral range from 350 to 1 000 nm,  giving a total  
close to 350 000 wavelength points per spectra. This huge number of 
thousands of stellar spectra, required to shape the table,
demands a lot of computer calculation time which damps the employment of
our techinque for the analysis of big databases,  such as 
Polarbase[http://polarbase.irap.omp.eu/]
<cit.>, where  hundreds of observed spectra of cool stars can be found.

To reduce the required number of stellar spectra to invert the MZSs, 
in this work we implemented 3 different regression methods based in machine 
learning algorithms.  All these three regression models belong to the 
supervised paradigm and are included in the Scikit-learn
software package <cit.> of the Python language. 
The employed regressions are: 1) The Bayesian Ridge Linear Regression (BR),
2) The Support Vector Machine (SVM), and 3) The Multi Layer Perceptron Artificial Neuronal 
Network (ANN). We next briefvly describe each one.









Regression analysis is a valuable mechanism to model and analyse data, constituting a form of 
predictive modelling technique used to explore the relationship between two (or more) variables: 
an independent one, the input (also known as feature or predictor), and a dependent one, the target.
Its goal is to  fit a curve to the data points in such a form that the sum of 
differences between the data points and the curve is minimal <cit.>. 
We next briefly describe each of the regressor models employed in this work.
 



 §.§ Bayesian Ridge


Linear Regression is one of the most frequently used modeling techniques:
It establishes  a linear relationship between the 
dependent variable  (y) and one or more independent variables (x⃗):


    y (x⃗,w⃗)=w_0 + w_1x_1 + w_2x_2 + ... +w_nx_n,

where the coefficients w_i are known as the parameters of the 
regression model. The first coefficient w_0 is the so-called 
interceptor or “bias” of the model. In our case, the input variable 
are the MZSs and the target can be any parameter
of either the stellar atmospheric model or the magnetic model 
(as for example T_ eff or H_ eff).

The linear  model that best fits the data is commonly calculated using the Least
Square Method and when this is the case, the performance of the
calculated model is measured by the R-square metric <cit.>. It is important to notice 
that this kind of regression require that a linear relationship between both
independent and dependent variables exists. It is also noteworthy the fact that
outliers terribly affect the regression line and, therefore, the estimated values.


If the power of the predictor on the estimated regression equation is bigger
than one, then the regression is known as polynomial. In that case, the regression
line is not a straight line but a curve that fits to the points. A higher power on
the independent variable allows a better fit for more complex datasets.

However, when the data experiences from multicollinearity, meaning that
it has multiple, highly correlated predictors, then the Ridge Regression technique
is used. When multicollinearity occurs, the Least Squares estimates are
unbiased but the variances are large, which in turn results on a higher error
on the predicted values. To correct this error, Ridge Regression adds a bias to
the regression estimates and it solves the multicollinearity problem through a
shrinkage parameter λ :


    min || y-x⃗w⃗ ||_2^2 + λ ||w⃗||_2^2


The shinkage parameter, λ≥ 0, penalizes the size of the coefficients on the 
original equation and it is the direct regulator of the amount of shrinkage, a process also
known as regularization <cit.>.

The regularization of the parameters, can be done using the Bayesian Regression techinque: it
introduce diffuse priors, also known as noninformative priors, which are 
probability distributions expressing general information. Using the Bayesian approach,
the regularization of the Ridge Regression model is equivalent to 
finding a maximum solution for a Gaussian prior (the previously mentioned diffuse prior)
with precision λ^-1.  However, instead of setting lambda manually, it is possible to 
consider it as a random variable to be estimated from the data if the predicted output y
is assumed to  follow a Gaussian distribution around x⃗w⃗ :


    p(y | x⃗, w⃗, α) = N (y | x⃗w⃗, α),

where α is a random variable and it needs to be estimated from the dataset.
A particular case of Bayesian Regression is the Bayessian Ridge Regression, where 
the prior for the parameter ω is given by a spherical Gaussian distribution, 
as in:


    p(ω | λ) = N (ω|0,λ^-1I⃗),

where I⃗ is the identity matrix. It is during the fit of the model
that the set of parameters ω, λ, α 
are estimated jointly <cit.>.



 §.§ Support Vector Machines


Support Vector Machines is a learning method used mostly to classify data and to detect
outliers, but  it can also be used for regression purposes. The SVM offers important 
characteristics in his performance, as for example high effectiveness in high 
dimensional spaces, usefulness even when the number of samples is lower than 
the number of dimensions and the memory efficency. 
In the SVM approach, data is plotted as points in a n-dimensional space
and the outcome of the algorithm will be an  optimal hyperplane that best 
categorizes the data clases  (for calssification purposes) or best predict 
the tendency of the data in order to predict future unkown data 
(regression model) <cit.>.


In the specialized literature, two main modalities are reported to perform regression 
tasks through the use of the SVM model, which are: the ϵ-support vector regression (ϵ-SVR) 
and the v-support vector regression (ν-SVR), <cit.>. 
In this work we have used the former approach ϵ-SVR:
Let X⃗ = (x⃗_⃗1⃗,x⃗_⃗2⃗,..,x⃗_⃗m⃗)
be a set of observations of the same independent variable, 
and let y⃗ = (y_1,y_2,...,y_m) the corresponding dependent variables. Also, let w⃗ be the weight vector, 
let b be a scalar, let Φ(x⃗) be a non-linear function, let ξ_i  be slack variables, and
let C>0 be their associated parameter. If ϵ>0 is the associated parameter with the 
ϵ-insensitive loss function, the ϵ-SVR solves the following optimization problem:


    min _w,b,ξ,ξ^*1/2w⃗^⃗T⃗w⃗ + 
     C ∑_i=1^mξ_i +  C ∑_i=1^mξ_i^*,

subject to

    w⃗^⃗T⃗Φ(x⃗_⃗i⃗) + b - y_i, ≤ϵ+ξ_i, 
    
    y_i -w⃗^⃗T⃗Φ(x⃗_⃗i⃗) -b, ≤ϵ+ξ_i^*, 
    ξ_i,ξ_i^* ≥ 0,i=1,...,m.


Similarly, the dual problem is defined according to α⃗_⃗i⃗:

    min _α, α^*1/2 (α⃗ -α⃗^*)^T -
     Q(α⃗ -α⃗^*)

subject to 

    e⃗^⃗T⃗ (α⃗ -α⃗^*) = 0 
    
    0 ≤α_i, α_i * ≤  C,      i=1,...,m

where the e^T is the identitity vector and Q_ij = K(x⃗_⃗i⃗,x⃗_⃗j⃗) ≡Φ(x⃗)^T Φ(x⃗)  
is the kernel function, 
which can be linear, polynomial, sigmoid or a radial basis function, among other options. 
The dual problem is solved by the Lagrange multipliers, and the resulting function is <cit.>:


    f(x⃗) = ∑_i=1^m (-α_i +-α_i^*) K(x⃗_⃗i⃗,x⃗_⃗j⃗)+b.






In our case, i.e. data analysis of MZSs,  we found after  testing differnet kernel fucntions 
that it was with the linear function that we  obtained the best predictions results.




 §.§ Artificial Neuronal Network



Artificial neural networks are inspired by physiological knowledge of the organization of the brain. 
They are structured as a set of interconnected identical units known as artificial neurons, and the 
interconnections are used to send signals from one neuron to the others, in either an enhanced or 
inhibited way. This enhancement or inhibition is obtained by adjusting connection weights. In an 
artificial neural network, a new observation causes an update of the network weights, which means 
that the network learns.

Let x_i ∈ℝ^n be the i-th observation from a set of m of them (corresponding to 
the same independent variable), and let x_ij∈ℝ be the j-th value of the i-th 
observation. From these m values and a previously established threshold value θ, it is 
possible to model a simple perceptron, which is an artificial (and very effective) model of 
a physiological neuron <cit.>. 
It is considered that the j-th value of the i-th observation 
(corresponding to a dendrite of a physiological neuron) has an associated weight w_ij∈ [0,1] 
when x_ij  becomes the entrance to the central part of the neuron (which corresponds to the soma 
of the physiological counterpart). The weight with value 1 is associated with the threshold value 
θ.

The process that occurs in the soma of the artificial neuron, is what gives the power to the 
machine learning algorithms based on artificial neural networks. This process involves a sum of the 
products and a function of activation of the neuron. Each product is obtained by multiplying the 
value of the neuron x_ij with its associated weight w_ij while the activation function
Θ is typically a sigmoid function or the tanh function. The result of this process, which 
is the output of the artificial neuron for of the i-th observation, is shown in the following expression:


    output(i)= Θ [  ( ∑_j=i^n w_ij x_ij) - θ ].








The simple perceptron is the fundamental building block of virtually all models of artificial neural 
networks of the state of the art (with some singular exceptions represented by unconventional models, 
such as morphological and Alpha-Beta models). The choice of the ways of interconnecting the simple 
perceptrons gives rise to different network topologies, among which stand out, undoubtedly the models
known as fully connected layered feedforward networks.


This type of models are organized in layers: an input layer, an output layer and one or more hidden 
layers. One of its important features is no connection or feedback between neurons of the same layer,
but the connections, as the name implies, are all feedforward; that is, a neuron located in a certain
layer, only connects with one or several neurons of the adjacent layer forward <cit.>.

The Multi Layer Perceptron (MLP) is a fully connected layered feedforward network that is trained 
with the method known as backpropagation, whose detailed description can be consulted in 
<cit.>. The MLP can be used as a classifier or as a regressor, and in this article 
we use a regressor version of the MLP known as MLPRegressor, which is trained with backpropagation 
using identitity as with the identity function as the activation function in the output layer, 
that corresponds to the output y of our application <cit.> .



 §.§ Regression model parameters

In practice, the implementation of each regressor model is defined by a 
set of hyperparameters, as for example the activation function  or the value 
of the penalization coefficients, etc. The choice of each hyperparameter 
for the three regressors used in this work are all included in Appendix A.

For the choice of the different hyperparameters in each regression model
we tested the inversion effectivity by varying the values of many  hyperparameters;   
the final values that were selected (listed in Appendix A), were those that offered  
better inversion performance for each regressor.
 
Finally, for the cases of the  SVM and the ANN regressors, the associated 
learning algorithms are sensitive to the data distribution. Both regressors assume 
that the data is centred around zero and is close to a standard normal distribution.
For this reason, it is requiered to apply a standardization of the data before the 
traning is performed. We have thus applied a pre-process to the data prior to the 
training process through the StandardScaler  utility included also in the 
Scikit-learn software package.




§ INVERTING MZSS

 
In this section we will employ the three different regression models
mentioned above to estimate the H_ eff using a supervised training. The reason
for using three different models is to determine to what extent
the results depend on the choice of the model, and thus to choose 
the most adequate one.





 §.§ Inverting noiseless MZSs


The MZSs that we will employ in this section are those previously established 
in , where we considered the case of a cool star adopting 
from the Atlas9 grid <cit.> an atmospheric model with 
T_ eff = 5 500 K, log g = 4.0, solar metallicity, zero macro-turbulence 
(v_ turb = 0.0) and a micro-turbulent velocity of ξ = 2.0 km s^-1.
Besides, we considered a slow rotator case fixing vsini = 10  km  s^-1. 
For the synthesis of the Stokes profiles we employed the code Cossam 
<cit.>. All the atmospheric parameters were kept fixed, but 
we allowed to randomly vary the 7 parameters that describe the configuration 
the dipolar tilted eccentric model (see above). The inversions showed that none of the free 
parameters could be retrieved being only possible to infer H_ eff. The 
fact that none of the free parameters of the dipolar model could be retrieved 
is due to the existent degeneracy among different combinations of parameters 
that can produce the same observable: H_ eff (see Fig. 10 in ). 
For this reason, the machine learning algorithms will be trained to predict 
only one target, the H_ eff.  We point out that our goal is to 
measure the stellar longitudinal magnetid field from spectropolarimetric snapshot 
data type,  and not to map the strength of magnetic field over the surface.



As first step, we must determine the (small enough) number of MZSs that 
we will use to train the regression models, keeping
in mind that it is required to obtain satisfactory and acceptable predictions
(inversion results).  











For this purpose, we first used only 50  MZSs to train the algorithms 
and then we tested each regression model inverting a representative set of  1  500 
MZSs. We then repeated the exersise but increasing the number of MZSs used 
in the training process by steps of 50 until reach a maximum of 500. The performance 
for all the regression models is shown in Fig. <ref>. 
The three models showed similar results, obtaining very accurate determinations
of the H_ eff even  when we considered a very small number of  MZSs in the 
training process: using 100 MZSs the obtained RMSE is inferior to 1 gauss. The general
tendency is that the errors decrease smoothly as we increase the number of  MZSs
employed in the training.  Based in these results, we have decided that training
the algorithms with 300  MZSs is sufficient for a proper determination 
of the stellar longitudinal magnetic field. Thus, hereafter all the inversions tests
will be performed using 300 MZSs in the training process. In  Fig. <ref> 
we show the individual inversions of the 1 500 MZSs for the three regression models. 


The excellent results of Fig. <ref> show the great
utility of the employment of machine learning algorithms in the inference of 
H_ eff from the MZSs. As comparison, in  we constructed  
7 500 spectra models and the RMSE of those inversions was 3.95 G, while with 
the best regression model (BR) we now obtain a RMSE of 0.41 G requiring only 
300 spectra models. In other words, using this regression model the inversions accuracy
is increased by one order of magnitud, and at same time the required time
to synthetize the stellar spectra models  (number of MZSs) is reduced by a factor of 25.
It is important to remark that the requiered time for the training of the algorithms is very fast :
using a single procesor the cpu-time  requiered was of 1, 5 and 120 seconds for the
BR, SVM and ANN regressors respectively; moreover,  in the case of BR and the ANN 
the time can be reduced since the training process be parallelized.

Unfortunately, the effectivity of the regression models is considerably reduced when 
noise is taken to account. To illustrate this, we repeat the same exersise
of Fig. <ref> but this time including 5 different noise levels 
in the MZSs. The noise level,  see Fig. <ref>, 
is quantified as the standard deviation of the noise divided by the standard deviation of the noise-free MZSs.
In Fig. <ref> we 
show the inversion results only for the BR model, which was the best 
regressor for the case of ideal MZSs without noise. (We note that our
definition of noise level is consistent with the work of <cit.>, so we can 
latter on make a proper comparison of the results.)

From the results of Fig. <ref>, it is clear that for the inversion 
of noise-affected MZSs,  the case corresponding to real data, the predictions are 
no longer good enough, even for the lowest noise level.  For these
tests, the training process was performed using noise-affected MZSs  
(with their respective noise levels), because in the contrary case 
– training with noiseless MZSs–, the errors were even slightly higher.
Besides, it is also worth noticing that we have verified that
the number of MZSs to include in the training process does not
depend on the noise levels, i.e., to consider 300 MZSs in the 
training process is good enough for all noise levels. 









Based on the obtained results when noise is included, 
we propose a two steps process to improve the efficiency of the
MZSs inversions. These two steps consist in a MZSs noise-cleaning and in an iterative 
inversion methodology.







 §.§ Denoising profiles with PCA




To facilitate the analysis of observed data, it is desirable to disposse of some
noise treatment in order to give more credibility to the results.
In the case of stellar spectropolarimetric data, the establishment
of the multi-line (LSD or MZS) profiles is already a pre-processing step 
of the  data analysis. However, the multi-line profiles are also
affected with noise, which in turn degradate the accuracy of the 
inversions. Thus, we will subsequently refer as noise treatment 
to any process applied to the multi-line profiles in order to deal with 
the noise problem.


Noise treatment, in this sense, can be implemented using for example 
any of two following approaches. On the one hand, the sparsity representation 
of data samples has recently begun to  be implemented in the analysis of 
solar and stellar spectropolarimetric data with very good results 
<cit.>. On the other hand, the principal 
components analysis (PCA), first introduced by <cit.> for the study of solar magnetic
fields, is a more commonly employed technique in the  analysis of both, solar and stellar 
spectropolarimetric data. Both approaches implement a sort of space dimensionality reduction 
allowing to diminish the impact of the noise in the retrieval of the model parameters.






We will employ PCA  to perform a “noise-cleaning” of the MZSs: In the previous section we decided 
to use 300 MZSs for the training of the regression algorithms. Applying the  
Singular Value Decomposition to this same set of MZSs we obtain in turn 300 eigenvectors. It is 
well known that the eigenvectors are a matemathical base of the original set of MZSs, such 
that any MZSs can be obtained as a linear combination of those eigenvectors <cit.>.
Nevertheless, it is also possible to consider not the entire set of eigenvectors but
only the first ones, because in fact the first eigenvectors are those that contain 
most of the usefull information about the original MZSs  <cit.>. 
Accordingly, we considered only the first ten eigenvectors to perform a noise-cleaning:




    MZS_cleaned≡∑ _i=1 ^10α_i    ev_i,

where α_i denote the scalar coefficients and ev_i  the eigenvectors.
The coefficients can be obtained through a dot product between each eigenvector 
and the MZS.




In Fig. <ref>, we illustrate the  cleaning proces of  MZS profiles considering the 
same 5 levels of noise previously used (from 0.2 to 1.0 in steps of 0.2). The solid lines in 
the upper panels correspond the noise added profiles, while in the lower panels they 
represent the respective cleaned profiles. As comparison, the dashed lines in both, 
upper and lower panels, correspond to the original noise-free profiles. It can be noticed
that the noise cleaning process works well even for the highest noise level. 
However, as expected, the higher the noise level, the greater the difference
between the original and the cleaned MZSs are found.


We applied the described process to the full set of 1 500 MZSs
and repeated the same inversion test of Fig. <ref>,
but this time implementing a noise cleaning prior to the inversion of the profiles. 
Surprisingly, the results are not as good as expected, being only 
slightly better to the case without cleaning process, 
as shown in Table <ref>. For this reason, we propose a complementary 
step in the data analysis  to improve  the inversion
efficiency.










 §.§ Iterative inversion process







Let us consider only one given MZS from which we would like to infer H_ eff, 
and let us also consider a given fixed noise level, for example 0.4. 
In the previous section we showed how to perform a noise cleaning of the 
profiles. We consider now a scenario in which we could repeat this process several times.
For each iteration, the random noise will affect differently the MZS and subsequently
the noise cleaned profiles will all be  slightly differents. In Fig. <ref>
we illustrate how the cleaned profiles at each iteration differ 
from the original noise-free profile. For example, the maximum amplitude for the 
bottom left cleaned profile is higher than the original one, contrary to 
the botom central cleaned profile where the maximum amplitude is inferior to the 
original one. These differences, even if very small, are enough to produce a different 
result when infering the longitudinal magnetic field. 
Nevertheless, we know that all the cleaned profiles have the same origin in the 
sense that all were established from the same noiseless data (dashed lines in Fig. <ref>). 
It is thus expected that the inversions derived from these cleaned profiles will be close to the real value of
H_ eff, and in fact they are, showing a distribution centred around
the original value of H_ eff.

















To illustrate our proposed methodology,  in Fig. <ref> we show two examples
of iterative inversions. We have considered two MZSs, one with an  H_ eff
of 238.2 G (left panels) and one with a much weaker magnetic  
longitudinal field of 5.2 G (right panels). From top to bottom, the three panels show 
how the inversion distributions varies with the number of iterations: 
the larger the number of iterations performed, the more the histograms will 
approximate a normal Gaussian distribution (solid line in black). The centroid of the fitted  
Gaussian distribution, indicated with a vertical line, corresponds to the inferred 
value of H_ eff in each histogram.

In order to inspect more carefully the dependence of the results on the number
of iterations, we tested the proposed inversion process using the full set of 1 500 MZSs.
Additionally, it is interesting to compare the performance of the three regression models 
introduced above. We recall that we are considering the same noise level of 0.4
for all the iterations and for all the MZSs. 
The RMSE of the inversions for the full set of profiles as function 
of the number of iterations is shown in Fig. <ref>. 
All the three regressor models show very similar behaivors, 
with a considerable improvement of the results specially for the first hundreds iterations:
When no iteration is applied the RMSE are 24.6, 29.7 and 26.3 G for 
the BR, SVM, and ANN respectively, and after 100 iterations the RMSE
is reduced to 16 G for the three models, and  at 1000 iterations the RMSE is of only 6 G. 
After around one thousand iterations, the improvement of the inversions becomes restrained.
With the present results,  we consider that apply 3000 iterations is good enough for the 
addopted inversion strategy.























Now, concerning the incertitudes associated to our measurements of H_ eff,  in  principle,
it is possible obtain an error estimation based in the the histograms of Fig. <ref>. 
Nevertheless, by doing it we would be highly 
overestimating the errors: The standard deviaton of the fitted Gaussian distribution for 
the bottom left panel is 160 G, but in fact the error is much more smaller, being of only 4.1 G. 
Simarly, an overestimation of the error is also obtained in other example of Fig. <ref>, 
at bottom rigth panel, where the standard deviation of the fitted distribution is 13.8 G but the 
real error is of only 1.2 G. Thus, in order to obtain realistic estimations of the errors, 
we next procced to inspect  the inversions over the full sample of 1 500 MZS.



 §.§ Caracterizing the technique









We will now show the performance of our proposed method considering differents 
noise levels; we tested the inversions using the same set of 1  500 MZSs as before,
and in all cases we applied 3  000 iterations. 
To quantify the accuracy in the results we have employed
the Mean Absolute Percentual Error (MAPE), defined as:





    MAPE=100 * (   |  H_ eff^original -   H_ eff^inferred  | /|  H_ eff^original |).



Additionaly, in order to underline the importance of including a noise treatment in the 
analysis of spectropolarimetric data, which in  our case is done through the noise 
cleaning of the MZSs followed by the iterative inversion procedure, we also included 
the results when the profiles were inverted without any noise treatment, 
which we labelled as single inversions.






In Fig. <ref>, we show the variation of MAPE
as function of the noise level. The three regressor models
(BR, SVM, ANNN) have the same performance as function of the noise level:
When no noise is included the MAPE is extremely low (0.1%), 
but the inclusion of even a very small noise level of 0.1  lead 
to a (relative) considerable increament of the MAPE (4%). 
However, the precision of the inversions for the  models 
will then only decrease moderately to finish with a MAPE of 11% for the highest 
noise level of 1. Finally, the case of traditional inversions –noted simple 
by the line with square symbols–, is the most affected one, varying from 13% for the 
lowest noise level of 0.1 up to  47%  for the highest noise level of 1. 



One conclusion from this test is that with our proposed
methodology -PCA noise cleaning combined with iterative
inversions-, the three regression models are robust methods to
determine H_ eff from MZSs that are affected by the noise. 
Another important conclusion is that the single inversions
are highly affected by noise, being the respective errors 
much higher than those obtained with any of the regression models.
This is particularly evident for the case of highest noise level 
of unity, where the MAPE of the single inversions can reach 
almost 50%, while for all three regressors, for this same noise level, 
the MAPE is close to 10%.

Similar results have been previously presented in  <cit.>.
In this work the noise treatment is done adopting a sparsity
representation of the data using an Ortogonal Match Pursuit (OMP) algorithm.
In Fig. 10 of that article, the authors presented an analogous test to the 
one just discussed here, comparing the inversions of noise affected profiles
using the OMP approach with the ones obtained using the traditional 
centre-of-gravity  (COG) method. The latter, is nowadays
the usual method employed when measuring  H_ eff from real data and it
is relevant in the sense that it does not consider any noise treatment. 
The authors have found that when considering a noise level of 1,
using the COG method the MAPE is close to 50% (the same as in our single case), 
while using the OMP approach the MAPE decreases to less than 20%.  
Both results, theirs and ours, stand out the importance of implementing a noise 
treatment in the inversion of polarised multi-line profiles.


The iterative process described in this section is a time-expensive step
in our approach since it takes 3.75 hrs using  8 processors in parallel. Nevertheless
this relative long time is requiered in order to achieve the precision errors showed
in Fig.  <ref>, and of course this time can be reduced if more processors
are at disposse. 




§ APPLYING THE METHOD TO REAL DATA





 §.§ HD 190771

In this section we will apply our method to real data. We have obtained from 
the public database PolarBase, 
spectropolarimetric data of two cool solar like stars, with
very similar physical properties to the Sun.
The first one is HD 190771, whose atmospheric parameters are:  T_ eff = 5834 ± 50 K, 
log g = 4.47 ± 0.03, M = 0.96 ± 0.13   M_, M/H=0.14 and 
vsini  =  4.3   km s^-1 <cit.>. 
For the synthesis of the stellar spectra, we took as starting point the closest model of 
the grid of Atlas9 atmospheric models (T_ eff = 5750 K, log g = 4.5 and  M/H=0.0),
to then extrapolate to the exact values of T_ eff, log g  and 
metallicity <cit.>.

The synthesis of the spectra covered a wavelength range of 369-1010 nm in steps
of 1 km s^-1. 
As mentioned in , when constructing the MZSs 
our only criteria for the line selection  is a minimum ratio of line-depth to continuum, 
which for this work we have fixed > 0.1, resulting in a total number of individual
lines of 23 172.
The established MZS for the intensity, circular polarisation 
and null profiles are shown Fig. <ref>. 















The next step is to fit the observed broadening for the intensity Stokes profile. 
For this propose, we set as free parameters the projected rotational velocity of 
the star (vsini) and the microturbulence velocity  (ξ):  
We calculated a sample of 50 stellar spectra varying vsini=[0.5,10] and ξ=[0,5]
km s^-1. We then  trained an algorithm to predict simoultaneously 
both parameters, obtaining as  best fit values vsini=5.97 and ξ=2.49 km s^-1.
In Fig. <ref> we show the obtained fit of 
synthetic profile to the observed one. 
In the adjustement of the intensity  MZS profile, we are considering that all the other 
broadening processes such as the magnetic, instrumental or macroturbulent ones, can 
be well represented by the interplay between only these two free parameters.
For this reason, the values of
vsini  and ξ derived from this fitting are not real measurements of these
two physical parameters.




In summary, we use an optimal atmospheric model (T_ eff=5834 K, 
log g=4.47, M/H=0.14, vsini=5.97 km s^-1 and 
ξ=2.49 km s^-1)  to synthesize a set of 300
MZS that we will use to train the regression algorithms and to subsequently

determine H_ eff for this star.






In order to follow a consistent MZSs inversion procedure
with our previous tests, it is required to determine the standar deviation of both, 
the noise-free MZS (σ_MZS) and the noise (σ_noi).
The latter was determined using the edges of the circular MZS profile 
–indicated by the dashed vertical lines in Fig. <ref>–, 
because in these regions there is only noise.
For the former, we considered only the central region 
where is present the polarised signature. 
We then applied  a noise cleaning process to the polarised MZS as described in 
the precedent section. Once we have obtained the cleaned MZS, we calculated
the standard deviation from this noise-free profile. Finally, the noise level is
determined as before, NL = σ_noi  / σ_MZS, obtaining a 
value of 0.15.

In order to apply the iterative inversions procedure, in each of the 39 
echelle orders we added random white noise to the observed circular polarised 
spectra. The amplitude (A_i) of the added noise in each order is different,
and it is given by :




    A_i = σ_i^ v / α,

where σ_i^ v is the standard deviation of the circular polarised spectra 
in  the i-order, and α is a constant factor for all orders. It is through the
value of α that the NL can be controlled in the MZS after the 
addition of random white noise. In Fig. <ref>
we show how the noise level of the MZSs varies as function of the α value
for this particular observation of HD 190771: 
For very low values of  α –corresponding to high amplitudes in the added noise–,  
the NL is very high,  but as α increases the NL decays very fast, to finally 
become stable  at around  NL ∼ 0.15, which correspond to the noise level 
of the MZS when no noise was added. 
In this plateau phase of the NL, corresponding to α values from 0.3 to 2, 
we are interested in the  minimum (α = 0.3), because the lower is the value of α 
the higher will be the added noise, i.e, the higher will be the variations in the 
iterative establishment of MZS profiles. 
In Fig. <ref>, we show some examples of the of MZS profiles 
for the HD 190771 data considerig α=0.3.














Finally,  we considered 3  000 iterations to produce the histogram of the inversions results.
In Fig. <ref> we show the distribution of the inversions,  as well as the fit 
of the normal Gaussian distribution, from which we derived a value of  
H_ eff = -6.3 ± 0.3 G. Using the same dataset, <cit.>
have reported a value of H_ eff = -9.8 ± 0.3 G. 



The value of the error in our measurement of H_ eff 
is based on the MAPE results of Fig. <ref>. 
Considering that in our inversions we employed the BR regressor and
given that the  noise level of the data is 0.15,   
it corresponds to a MAPE close to 5%, i.e., an error of ± 0.3 G. 
It is worth noticing  that the same 
error is reported by  <cit.> despite the fact that in the inversion of the LSD 
profile, the authors they do not implement any noise treatment. As a reference,
for the same NL of 0.15, the MAPE associated to single case 
(inversions without noise treatment) is 15%, i.e., 
the error should be around ± 1.5 G. It seems thus  that the errors
reported by classical techniques that do not consider noise traetment could be 
seriously underestimated. 


Now, if we consider the estimation that we found of 1.5 G for the error  in the 
LSD measurement of H_ eff,  then it turns out that the value of  H_ eff reported by 
<cit.> and the one we found,  would be in agreement at a difference of
2σ. On the contrary, if we consider the error of 0.3 G reported by the authors, then 
the difference could not be explained by the incertitudes associated to the measurements, 
and in this case the discrepance in the results deserve
a much deeper analysis, which is beyond the scope of the present study.





 §.§ HD 9472






The second example is the star HD 9472, whose atmospheric parameters are
T_ eff=5867 ± 44 K, log g=4.67 ± 0.02, 
M=1.29 ± 0.19   M_, M/H=0.0 and  vsini=2.2   km s^-1 <cit.>.


We have followed the same procedure as in the previous case, taking as point of departure
the closest atmospheric model from the Atlas9 grid, namely, T_ eff = 5750 K, 
log g = 4.5, and  M/H=0.0, to then extrapolate to the exact 
values of T_ eff and  log g. Adopting the same threshold of the line-depth to continuum as before (> 0.1), 
the total number of individual lines is 23 895 for this star.  We then established the MZSs,
shown in Fig. <ref>.

The next step was to reproduce the observed broadening in the intensity Stokes MZSs.
We obtained as best fit values vsini=5.08 and ξ=1.81 km s^-1,
which produced, as in the precedent case, a very good fit between the synthetic and
the observed profiles, as showed in  Fig. <ref>.









The NL was then calculated  for the circular MZS, obtaining a
value of 0.19.  We subsequently inspected the variation of the noise level as function
of the amplitude of the random noise added to the echelle orders, 
allowing to determine that setting α =0.4 will correspond
to MZSs with the same NL as the original data of 0.19, see Fig. <ref>. 
Finally, we  iteratively stablished 3  000 MZSs 
to consequently fit a Gaussian distribution to the histogram of the inversions,  
shown in Fig. <ref>. The centroid of the fitted distribution 
corresponds to a value of H_ eff = -3.3 ± 0.2 G, which in fact 
is almost the same that the one reported  by  <cit.> using the same data:
H_ eff = -3.5 ± 0.4 G. Once more, it seems that the error
reported by the authors is underestimated: based in Fig. <ref>,
for this noise level the error should be around 17%, i.e., 0.6 G.




§ CONCLUSIONS



We have presented a new study related to the measurement of stellar 
longitudinal magnetic fields, from high resolution spectropolarimetric data,
analysing multi-line profiles. 
Our main goal was to develope a general method which is constrained nor by 
the regime of validity of the weak field approximation
neither by the line autosimilarity normally assumed in the classical 
methods employed nowadays in the measurement of H_ eff. 

Our technique is thus based on a theoretical radiative transfer approach where  
we produced a synthetic set of Stokes profiles, each one corresponding to a different
configuration of the magnetic fields over the stellar surface. Considering the best 
possible scenario of noise-free spectra, we have showed that the 
use of machine learning algorithms (MLA) is a key step to reduce the number of 
synthetic profiles required for an acceptable accuracy in the correct 
determination of H_ eff: With the use of MLA and 
considering only 50 stellar spectra, we obtained a better value of the RMSE for the inversion
results that when we considered 7  500 spectra but the inversion strategy was based 
in a look-up table, i.e.,  no MLA was applied <cit.>.

We also showed that the results are considerably degradated when noise is included; 
nevertheless this is not restricted to the MLA but in general to any method dealing with 
the inversion of noise-affected profiles or signals. We thus proposed a data noise pre-processing
consisting of two steps: 1) To use PCA to perform a noise-cleaning of the multi-line  profiles 
and 2) An iterative inversion procedure.  Applying this data analysis process, we achieved a considerable
improvement in the inversion accuracy, confirming that it is desirable  to include a 
noise treatment  in the analysis of multi-line profiles in order 
to get more confidence in the derived values of H_ eff. Very similar results
about the  impact of the noise  in the inversion of multi-lines profiles were independently found 
by <cit.>.


We implemented our inversion technique to real observations of two stars. In order to
similarly apply the MLAs to the observed data set (i.e. as closely as possible
as in the case of synthetic spectra),  
it was first required to reproduce 
the observed broadening in the intensity profile. We allowed the variation of
two physical parameters for this propose, namely, vsini and ξ, and 
we obtained very good fits, as shown in Figs. <ref>
 and <ref>. Of course, the  values of these two parameters
lack of any physical interpretation because we are assuming that all other
broadening mechanisms (e.g. magnetic, instrumental or any other) can be well
represented by considerig only vsini and ξ as adjustable parameters.

Thus, the atmospheric model employed for the analysis of each star consisted in  
the exact values of T_ eff, log g and metalicity –reported 
by other authors–, in addition to the values of vsini and ξ. 
Once these five atmospheric parameters were defined, we synthesized a set 
of 300 polarised stellar spectra that in turn were used to train the MLA
(the used training data sets can be found in <www.astrosen.unam.mx/j̃ulio/ML_mzs>).
One final step was to reproduce the iterative inversions methodology implemented in the
synthetic tests. For this purpose, we added random white noise in each 
order of the echelle spectrum. The amplitude of the added noise in each order 
is different and is controlled by one free paramter, labelled α in Eq. (<ref>). 
It is through the value of this parameter that it was possible to keep a constant noise level 
in  the MZSs established for the iterative inversions. Finally, considering the same 
number of iterations in the case of real 
data than in the tests,  we fitted  a Gaussian normal function to the distribution
of the inversion results. The centroid of the Gaussian distribution  determines the value of 
H_ eff for each data set (Figs. <ref> and <ref>).


We have used two stars as test cases to illustrate the direct application of 
the proposed methodology to real 
observations, and in this sense,  it is not the scope of this work
to compare the measured values of H_ eff with our technique
with the ones obtained by other methods. We will proceed to the comparison 
of the values of  H_ eff obtained with classical methods versus ours 
in a forthcoming article.
However, what we can highlight the fact that the uncertainties
reported in the measurement of H_ eff form the LSD profiles are most likely
underestimated. 
<cit.> have recently presented a method to establish
and analyse LSD profiles under a Bayessian framework. One of the
advantages of this approach is that it can estimate the intervals of credibility 
at each velocity point of the LSD profiles. Unfortunately, 
although the authors applied their method to three different stars, 
they did not report the respective measurements of H_ eff, preventing
the comparaison of their estimation of errors with other techniques.
As far as we know, no other work has addressed a detailed study 
on the expected uncertainties associated to the measurements of H_ eff 
from noise-affected LSD profiles.

The most relevant aspect of the work presented here is not that 
we have included estimations of the uncertainties in the measurements
of H_ eff, but the fact that we have introduced a new 
methodology for the analysis of polarised multi-line profiles
–including a radiative transfer theoretical framework in the synthesis
of the Stokes profiles–, 
using a strategy based on machine learning algorithms. 
We  expect  to apply our data analysis method to stars other 
than solar-type stars in incoming studies.





§ ACKNOWLEDGEMENTS

We are very greatful to the whole group of Laboratorio de Cómputo Inteligente del CIC-IPN
[http://www.cic.ipn.mx/sitioCIC/index.php/pre-rncnc] 
for very constructive and illustrative discussions
about machine learning methods. We would also like to express our acknowledgement to
the anonymous referee for all the comments that helped to improve the clarity of
this work.
JRV acknowledges support grants from CONACyT 240441 and
Supercomputo – LANCAD-UNAM-DGTIC-326.
Thanks go to AdaCore for providing the GNAT GPL Edition of its Ada compiler.



aa




§ HYPERPARAMETERS FOR EACH REGRESSOR

For the definition of each of the hyperparameters listed bellow, we refer the reader 
to Scikit-learn user guide. For consistence, we labeled each  parameter 
following the mentioned guide 
[http://scikit-learn.org/stable/_downloads/scikit-learn-docs.pdf].


 §.§ Bayessian Ridge

BayesianRidge(alpha_1=1e-06, alpha_2=1e-06, compute_score=False, copy_X=True,
       fit_intercept=True, lambda_1=1e-06, lambda_2=1e-06, n_iter=300,
       normalize=False, tol=0.001, verbose=False)



 §.§ Artificial Neuronal Network

MLPRegressor(activation='identity', alpha=0.05, batch_size='auto', beta_1=0.9,
       beta_2=0.999, early_stopping=False, epsilon=1e-08,
       hidden_layer_sizes=(100,), learning_rate='constant',
       learning_rate_init=0.001, max_iter=10000, momentum=0.9,
       nesterovs_momentum=True, power_t=0.5, random_state=25, shuffle=True,
       solver='lbfgs', tol=1e-08, validation_fraction=0.1, verbose=False,
       warm_start=False)
       


 §.§ Support Vector Machine

SVR(C=50.0, cache_size=200, coef0=0.0, degree=3, epsilon=0.1, gamma='auto',
  kernel='linear', max_iter=-1, shrinking=True, tol=0.001, verbose=False)

       

