Machine Learning Cosmic Expansion History
Deng Wang Wei Zhang
astro-ph.CO astro-ph.HE gr-qc



cstar@mail.nankai.edu.cn
Theoretical Physics Division, Chern Institute of Mathematics, Nankai University,
Tianjin 300071, China

cosmoszhang@mail.nankai.edu.cn
Department of Physics, Nankai University, Tianjin 300071, China

We use the machine learning techniques, for the first time, to study the background evolution of the universe in light of 30 cosmic chronometers. From 7 machine learning algorithms, using the principle of mean squared error minimization on testing set, we find that Bayesian ridge regression is the optimal method to extract the information from cosmic chronometers. By use of  a power-law polynomial expansion, we obtain the first Hubble constant estimation H_0=65.95^+6.98_-6.36 km s^-1 Mpc^-1 from machine learning. From the view of machine learning, we may rule out a large number of cosmological models, the number of physical parameters of which containing H_0 is larger than 3. Very importantly and interestingly, we find that the parameter spaces of 3 specific cosmological models can all be clearly compressed by considering both their explanation and generalization abilities.












Machine Learning Cosmic Expansion History
    Wei Zhang
     December 30, 2023
=========================================



Introduction.— During the past two decades, with gradually mounting data, a large number of cosmological observations such as Type Ia supernova (SNIa) <cit.>, cosmic microwave background (CMB) anisotropies <cit.> and baryonic acoustic oscillations (BAO) <cit.> have indicated that the universe is undergoing a phase of accelerated expansion, which is dubbed as dark energy (DE) by cosmologists. Recently, this mysterious phenomenon of cosmic acceleration has been confirmed, once again, by the galaxy clustering and weak gravitational lensing data from the first year (Y1) release of the Dark Energy Survey (DES) <cit.>. In general, to characterize the properties of DE based on observational data, theoreticians propose statistically two methods, i.e., parametric and non-parametric methods. For the former case, one can specify a DE model based on some phenomenological considerations or physical mechanisms, and then constrain its parameter space using various observations. For example, up to now, the simplest DE model is the cosmological constant (Λ) plus cold dark matter (CDM) scenario, namely the so-called ΛCDM model, which has phenomenologically one typical parameter (matter density ratio Ω_m) if we only consider the late-time evolution of the universe <cit.>. What one should do is to constrain ΛCDM using data and calculate the 1-dimensional confidence interval of Ω_m. The goals of constraints are: (i) for a given DE model, to determine the optimal values and statistical confidence intervals of key physical parameters; (ii) to rule out a class of DE models, which are not supported by statistical information criteria. For the latter case, without any particular parameterized form, one can reconstruct the evolution of DE or study the values of specific cosmological quantities starting directly from observations. Note that these methods must be more statistically complex than parametric ones. In the literature, there are several popular non-parametric methods such as principle component analysis (PCA) <cit.>, local weighted regression (LWR) <cit.> and Gaussian process (GP) <cit.>. PCA is widely used for reducing dimension and denoising of large datasets and it has been introduced to reconstruct the equation of state (EoS) of DE <cit.>. LWR, which is a k-nearest-neighbor-based method, is a general framework of least square regression. It has also been used for studying the dynamics of DE <cit.>. GP is recently used with high frequencies in the fields of astrophysics and cosmology <cit.>.
It is a stochastic process and its distribution is a distribution over functions with a continuous domain, e.g., time or space. It has been used for studying the evolution of EoS of DE <cit.> and reconstruct the CMB angular power spectrum <cit.>.

It is worth noting that previous researches all focus on testing the explanation ability (EA) of a given DE model in the usage of parametric methods, i.e., implementing the standard χ^2 statistics to estimate the parameter space of typical model parameters in light of all data <cit.>. Meanwhile, the same consequences occur in the usage of non-parametric methods, i.e., using all data, previous works all concentrate on reconstructing the underlying models with the best EA, which govern the dynamical evolution of cosmological quantities <cit.>.

However, for a given set of data, the generalization ability (GA) of a model learnt from some statistical algorithm is also very important from a point of view of machine learning (ML), which is the core of artificial intelligence and data science. Usually, in practice, a ML specialist only considers the GA of a model, but one statistician just studies the EA of a model. With gradually accumulating cosmic data and improvements of computing abilities, we think that cosmologists faces a challenge, i.e., developing and using some advanced statistical method to find a DE model, which has both good EA and GA for a group of datasets. Nonetheless, whether the EA and GA of a DE model are good or not is very hard to judge and quantify simultaneously. Furthermore, the fitted DE model from data may have no reasonable physical form. To overcome naturally both difficulties, we transform the above challenge into its dual question: fixing a specific DE model such as ΛCDM, how can we find out the overlapped region of two parameter spaces from its best EA and GA in light of ML techniques ? Therefore, in this work, our task is to obtain the overlapped and reduced parameter space based on traditional χ^2 statistics and ML algorithms for a given DE model and a given set of data.

Data.— As a performance, we adopt 30 cosmic chronometers (CC) lying in the redshift range (0, 2) as our data sample <cit.>, which are obtained using the most massive
and passively evolving galaxies based on the “ galaxy differential age ” method, and consequently are independent of any DE model.

Methodology.— Firstly, we perform χ^2 statistics to constrain ΛCDM using CC, and let the obtained parameter space via the Markov Chain Monte Carlo (MCMC) technique correspond to the best EA of ΛCDM. Secondly, through separating the whole sample into training set and testing set, we take seven different ML algorithms to learn CC data, and use the performance measure, mean squared error (MSE), to choose the optimal ML algorithm, which will be used for obtaining the parameter space corresponding to the best GA of ΛCDM via MCMC method. Furthermore, we can find out the overlapped parameter space for ΛCDM. Then, we shall repeat this process for non-flat ΛCDM (oΛCDM) and ωCDM models, where ω is constant EoS of DE.

Specifically, 7 ML algorithms are, respectively, ordinary least square (OLS), least absolute shrinkage and selection operator (LASSO) <cit.>, ridge regression (RR) <cit.>, elastic net (EN) <cit.>, linear support vector regression (LSVR) <cit.>, neural network (NN) <cit.> and Bayesian ridge regression (BRR) <cit.>. To choose the optimal ML algorithm, at first, except for OLS, we select a specific value for the hyperparameter according to our previous experiences, and take randomly 3-fold cross validation to train data and test the obtained model for each algorithm. Subsequently, minimizing the loss function (LF) on training set, we can get a model with best-fitting values of parameters, which shall be used for calculating MSE on testing set. Then, we implement the process of hyperparameter tuning and choose the best value of hyperparameter, which makes MSE on testing set minimal. Repeating this program 3 times and averaging 3 MSEs, we can have the best MSE, which is regarded as our score for each algorithm. Finally, we select the optimal algorithm which has the least score among 7 algorithms.

When applying ML to CC, we first build a bridge between ML and χ^2 statistics by identifying the LF in ML with the likelihood function in χ^2. Since one in the field of ML often deals with data without errors but cosmological data such as CC has errors, we then add errors of CC, for the first time, into the LF by fitting CC with a power-law polynomial (PLP), where we add errors of Hubble parameters H(z) at each redshift z into each term of PLP. Note that, to extract utmostly the information hiding in CC, we have taken a feature adding technique via a PLP to express H(z), since there is only 1 feature for CC using ML terminology. Subsequently, we propose a conjecture:

“ The model learnt from data is equivalent to the simplest one constructed from physical considerations ”.

For the case of CC, it can be shown as “ H_ML(z)= H_ΛCDM(z) ”  (or H_oΛ CDM(z), H_ω CDM(z)):

    ω_0+ω_1z+...+ω_n_dz^n_d=H_0√(Ω_m(1+z)^3+1-Ω_m),

where H_0, n_d and ω_i (i∈ [0, n_d]) denote the Hubble constant, the number of degree and the coefficients of PLP, respectively. To determine the 2-dimensional parameter space (H_0 versus Ω_m ) of ΛCDM, we must first determine n_d of H_ML(z), which equals the number of terms n_t of H_ML(z) minus 1. Then, we Taylor expand the above equation on both sides around z=0 and can naturally reexpress the coefficients ω_i of PLP with H_0 and Ω_m. Furthermore, inserting the reexpressed ω_i into the LF with optimal hyperparameters, we minimize the LF on the whole sample and thereby obtain the parameter space from ML, which represents the best GA of ΛCDM.
( The basic formula used in this work and all the mathematical details of 7 ML algorithms can be found in the supplementary materials.)










Results.— For 6 parametric ML algorithms, using CC, we find that the best degree n_d=2 for a wide range n_d∈ [1, 30] and this result is very stable for 3-fold data separation, and that the zeroth-order term of each PLP is H_0. Meanwhile, in the upper panel of Fig. <ref>, the optimal algorithm which has the least MSE 1.022 among 6 algorithms is BRR. Since we take conventionally OLS to confront cosmological models with observations, OLS should be used as as an important reference. We find that OLS has the largest MSE 2.957 and performs the worst at n_d=2. It is noteworthy that MSE of LASSO obeys a different variation tendency in comparison with those of 5 left methods, i.e., it suddenly decreases at n_d=4. Interestingly, its value at n_d=2 is still less than that at n_d=4 by a difference 0.162. The MSE of LSVR 1.291 is a little larger than those of EN and RR, which are very similar at all degrees.

For completeness and a naive comparison to 6 parametric methods, we also introduce the non-parametric error-back-propagation NN algorithm, which consists of 1 input layer (1 neuron), 2 hidden layers (30 neurons per layer) and 1 output layer (1 neuron). Note that for the first time in cosmology, we have successfully adding data errors into NN, which is very important to maximize the information utilization. After careful hyperparameter tuning, we find that BRR is still the optimal algorithm (see the lower panel of Fig. <ref>). Hence, using CC, we choose BRR as our best ML method to derive the generalization parameter space (GPS) of ΛCDM.

The relationships between H(z) and z from 7 ML algorithms only using CC and ΛCDM (here we use H_0=66.93 km s^-1 Mpc^-1 <cit.> and Ω_m=0.308 <cit.>) constrained by Planck CMB angular power spectrum data are shown in Fig. <ref>. We find that the best-fitting models of 6 parametric algorithms are very close to each other as well as ΛCDM and that of our non-parametric NN exhibits an interesting straight-line-like behavior. Furthermore, in Fig. <ref>, we present 1σ (%68) and 2σ (%95) confidence bands from our best BRR method, and obtain the first Hubble constant estimation H_0=65.95^+6.98_-6.36 km s^-1 Mpc^-1 from BRR via a n_d=2 PLP by only considering the GA of the model learnt from data. Although this new estimation is consistent with the indirectly global derivation by Planck Collaboration <cit.> and directly local measurement by Riess et. al. <cit.> at the 1σ confidence level (CL), its optimal value prefers the Planck's derivation. Note that this estimation is not only very compatible with our previous result H_0=67.38±4.72 km s^-1 Mpc^-1 based on GP at the 1σ CL by using the same CC data <cit.>, but also has larger uncertainties than it does. In addition, one can also easily find that the best-fitting ΛCDM from Planck data is very consistent with the underlying model from BRR at the 1σ CL.

As described above, optimizing the reexpressed LF from BRR, we obtain, respectively, the best-fitting values and 1σ uncertainties of parameters corresponding to GAs of ΛCDM, ωCDM and oΛCDM models. Specifically, we have H_0=65.07^+4.35_-4.06 km s^-1 Mpc^-1
, Ω_m=0.292^+0.094_-0.080 for ΛCDM, H_0=65.89^+7.68_-13.35 km s^-1 Mpc^-1, ω=-0.751^+1.096_-0.304 for ωCDM, and H_0=65.89^+8.47_-7.61 km s^-1 Mpc^-1, Ω_k=0.064^+0.363_-0.247 for oΛCDM. Subsequently, we find that: (i) although the values of H_0 from GAs of 3 models are consistent with those from their EAs at the 1σ CL, the best-fitting H_0 values from GAs are all less than those from EAs; (ii) the values of Ω_m, ω and Ω_k for 3 models are, respectively, compatible with Planck's results at the 1σ CL <cit.>; (iii) for ωCDM, considering GA, ω is weakly anti-correlated with H_0, which is different the case of EA; (iv) all the values of parameters have large errors due mainly to small data sample size.

Plotting the explanation parameter space (EPS) from χ^2 statistics together with the GPS from BRR, we obtain the overlapped and reduced parameter spaces of ΛCDM, ωCDM and oΛCDM models, respectively (see Fig. <ref>). Furthermore, we specify these narrow and overlapped regions at the 2σ CL as follows: H_0∈ [60.52, 72.14], Ω_m∈ [0.206, 0.454] for ΛCDM, H_0∈ [56.28, 74.49], ω∈ [-1.206, -0.335] for ωCDM, H_0∈ [56.86, 75.97], Ω_k∈ [-0.255, 0.560] for oΛCDM. From the view of statistics, we find that oΛCDM from EA has a stronger GA than ΛCDM and ωCDM from EAs, since the area ratio R between the overlapped contour and the EA contour of oΛCDM is larger than 50% (R<50% for both ΛCDM and ωCDM). From Fig. <ref>, one can also find that the 1σ interval of H_0 of overlapped parameter space of ΛCDM is smaller than those of ωCDM and oΛCDM. This indicates that the situation that the EPS of ΛCDM should become larger than before by adding a physical parameter ω or Ω_k to it will also occur to its overlapped parameter space. In total, making full use of CC, we have verified that the parameter spaces of 3 cosmological models can be obviously compressed (or reduced) by considering both their EAs and GAs.


Discussions.— Using CC, our best degree n_d=2 for 6 algorithms may imply that, from the view of ML, we have ruled out a large class of models where the number of physical parameters n_p>3 including H_0. This is why we just choose the ΛCDM (n_p=2), ωCDM (n_p=3) and oΛCDM (n_p=3) models to derive the overlapped parameter spaces. Based on the consideration that the ultimate aim using data to constrain DE models is to reduce ceaselessly their EPSs, discard all the impossible models and give the more accurate values of parameters of correct models than before, to a large extent, we can reduce the current EPSs if we also consider the GPSs of DE models. From Fig. <ref>, we know that BRR in all 7 ML algorithms is the best algorithm to characterize the GA of 3 different cosmological models for CC data. But if one is dealing with other datasets such as SNIa, LSS and CMB, the best ML method may not be BRR. Although overlapped parameter spaces of the above 3 DE models are all found out by using CC, we cannot ensure that the EPS and GPS of any given DE model with 1<n_p⩽3 have a common region.

Conclusions.— We have successfully introduced the ML techniques into the field of cosmology by considering simultaneously the EA and GA of a specific DE model. For the first time, we add the information of data errors into the usage of 7 ML algorithms via a PLP expansion when studying the background evolution of the universe.
We make a connection between ML and traditional χ^2 statistics by use of a physical example, i.e., using BRR to find out the overlapped regions of EPSs and GPSs of ΛCDM, ωCDM and oΛCDM models, respectively. For the first time, using CC, we find that the parameter spaces of 3 DE models can be clearly reduced by considering both their EAs and GAs. Very interestingly, from the view of statistics, we find that oΛCDM has a stronger GA than ΛCDM and ωCDM by performing χ^2 fitting to CC data.

Prospects.— With fast growing astronomical data, continuously improved data-storage capabilities and computing abilities, there is no doubt that ML will play a very important role in the researches of astronomy, astrophysics and cosmology. We expect that ML methods can extract new and useful information about underlying physical laws from future huge amounts of data.

Acknowledgements.— We thank Jing-Ling Chen, Wu-Sheng Dai and Yu-Xiao Liu for useful communications.



99

1
A. G. Riess et al. [Supernova Search Team], Astron. J. 116, 1009 (1998).

2
S. Perlmutter et al. [Supernova Cosmology Project], Phys. Rev. Lett. 83, 670 (1999).

3
C. L. Bennett et al. [WMAP Collaboration], Astrophys. J. Suppl. Ser. 208, 20 (2013).

4
P. Ade et al. [Planck Collaboration], Astron. Astrophys. 571, A16 (2014).

5
D. H. Weinberg et al., Phys. Rep. 530, 87 (2013).

6
  T. M. C. Abbott et al. [DES Collaboration],
  
  arXiv:1708.01530.

7
  S. M. Carroll,
  
  Living Rev. Rel.  4, 1 (2001)

8
 K. Pearson, 
 Philosophical Magazine. 2 (11): 559-572 (1901).

9
H. Hotelling,  
 Journal of Educational Psychology, 24, 417-441, and 498-520 (1933).

10
 W. S. Cleveland, 
  74 (368): 829C836 (1979).

11
 W. S. Cleveland and Susan J. Devlin, 
 Journal of the American Statistical Association. 83 (403): 596C610 (1988).

12
 M. Seeger, 
 International Journal of Neural Systems. bf 14 (2): 69C104 (2004).

13
D. Hurterer and G. Starkman, Phys. Rev. Lett. 90, 031301 (2003).


14
A. Montiel et al., Phys. Rev. D 89, 043007 (2014).


15
T. Holsclaw et al., Phys. Rev. Lett. 105, 241302 (2010).


16
M. Seikel et al., J. Cosmol. Astropart. Phys. 06 (2012) 036.




17
  D. Wang and X. H. Meng,
  
  Phys. Rev. D 95, 023508 (2017).

18
  D. Wang and X. H. Meng,
  
  Astrophys. J.  843, 100 (2017).

19
  D. Wang and X. H. Meng,
  
  Phys. Dark Univ.  18, 30 (2017).



20
  A. Aghamousa, J. Hamann and A. Shafieloo,
  
  J. Cosmol. Astropart. Phys. 09 (2017) 031.

21
M. Moresco et al.,  J. Cosmol. Astropart. Phys. 05 (2016) 014.

22
R. Tibshirani, 
Journal of the Royal Statistical Society. Series B (methodological). Wiley. 58, (1): 267-88 (1996).

23
A. N. Tikhonov and V. Y. Arsenin, Solution of Ill-posed Problems. Washington: Winston & Sons (1977). ISBN 0-470-99124-0.

24
H. Zou and T. Hastie, 
Journal of the Royal Statistical Society, Series B: 301-320 (2005).

25
H. Drucker et al., “ Support Vector Regression Machines ”, in Advances in Neural Information Processing Systems 9, NIPS 1996, 155C161, MIT Press.

26
W. McCulloch and W. Pitts, 
Bulletin of Mathematical Biophysics. 5 (4): 115C133 (1943).

27
M. Goldstein and D. Wooff, Bayes Linear Statistics, Theory & Methods. Wiley (2007). ISBN 978-0-470-01562-9.

28
N. Aghanim et al. [Planck Collaboration], Astron. Astrophys. 596, A107 (2016).

29
A. G. Riess et al., Astrophys. J. 826, 56 (2016)

30
  D. Wang and X. H. Meng,
  
  Sci. China Phys. Mech. Astron.  60, 110411 (2017).





















































































































