Reconstruction of latetime cosmology using Principal Component Analysis

Cosmological parameters are now constrained to much better precision than before Aghanim:2018eyx . This has been facilitated with significant improvements in observational techniques and the accessibility of various observational data. Data in the last two decades has confirmed that the dynamics of the Universe is dominated by a negative pressure component, known as dark energy. Dark energy is understood to be the component driving the observed accelerated expansion of the present Universe. One of the prime endeavors of modern cosmological research is to reveal the fundamental identity of dark energy, its exact nature and evolution.

It is not clear from the present observations whether dark energy is a cosmological constant Carroll1992 ; Carroll:2000fy ; Turner:1998ex ; Padmanabhan:2002ji or a time-evolving entity Peebles:2002gy ; Copeland:2006wr . The dark energy can be described by equation of state parameter $w=-P_{de}/\rho_{de}$, where $\rho_{de}$ is the energy density and $P_{de}$ is its pressure contribution. The form of the equation of state parameter $w$ of dark energy depends on the theoretical scenario being considered. A constant value $w=-1$ corresponds to the $\Lambda$CDM (cosmological constant with cold dark matter) model, whereas in case of time-evolving dark energy, the equation of state parameter can have different values Carroll1992 ; Padmanabhan:2002ji ; Peebles:2002gy ; Weinberg1989 ; Coble1997 ; Caldwell1998 ; Sahni2000 ; Ellis2003 ; Linder2008 ; Frieman2008 ; Albrecht2006 ; licia2008 ; nesseris2020 . We have little theoretical insight into these models, except for the $\Lambda$CDM model, which has a strong theoretical motivation. However, the standard ($\Lambda$CDM) model faces the problem of fine-tuning Carroll1992 ; Padmanabhan:2002ji ; Weinberg1989 ; Coble1997 . The observed value of the cosmological constant is found to be many orders of magnitude smaller than the value calculated in quantum field theories. Various alternative models have been proposed, which are based on fluids, canonical and non-canonical scalar fields. These models have the fine-tuning problem of their own, for instance, scalar fields require potentials which are specifically tailored to match observations Ratra:1987rm ; Linder:2006sv ; Caldwell:2005tm ; Linder:2007wa ; Huterer:2006mv ; Zlatev:1998tr ; Copeland:1997et ; PhysRevD.66.021301 ; Singh:2019bfd ; Singh:2018izf ; Bagla:2003prd ; Tsujikawa:2013 ; Rajvanshi:2019wmw ; Chevallier:2000qy .

Since the approach to understanding dark energy is primarily phenomenological, it is necessary to determine the model parameters to constrain and rule out models that are not consistent with data. Likelihood analysis is the most commonly used technique in cosmological parameter estimation and model fitting Singh:2019bfd ; Singh:2018izf ; Jassal:2009ya ; Nesseris:2004prd ; Nesseris:2005prd ; Nesseris:2007 ; Archana:2018 ; Sangwan:2018zpz . It is based on Bayesian statistical inference, where the posterior probability distribution of a parameter is determined with a uniform or variant prior function and the likelihood function. A combination of different data-sets, with likelihood regions complementary to each other allows a very narrow range for cosmological parameters. Thus the combination of various data-sets tightens the constraints on the parameters. This parametric reconstruction Chevallier:2000qy ; Linder:2002et ; Jassal:2004ej ; Gong:2006gs ; Mukherjee:2016eqj ; 2018PhRvD..98h3501V ; DiValentino:2017zyq ; licia2020_1 ; licia2020_2 ; licia:2013 assumes a functional form of the equation of state parameter of dark energy. However, it induces the possibility of bias in the parameter values due to the assumption of the parametric form. One elegant way to constrain the equation of state parameter is by the differential age method of cosmic chronometer technique, which can be used to find out the evolution of the Universe without any inference from a particular cosmological model licia2008 ; Jimenez:2001gg ; Moresco:2018ApJ ; Valent:2018 .

An alternative method is to reconstruct the evolution of cosmological quantities in a non-parametric fashion. Various statistical techniques have been adopted for non-parametric reconstruction of cosmological quantities Montiel2014 ; Gonz2016 ; Taylor2019 ; Porqueres2017 ; Diaz2019 ; Arjona2019prd . One of the promising techniques in non-parametric approaches is the Gaussian Process(GP) Valent:2018 ; Sahlen:2005zw ; Holsclaw2010am ; Shafieloo:2012ht ; Francesca:2019jcap ; Rasmussen:2005 ; Bonilla2020 , where we can create a multivariate Gaussian function with a determined mean and corresponding covariance function from any finite number of collections of random variables.

To accomplish the reconstruction of cosmological quantities, we use the Principal Component Analysis (PCA). A comparative study of different model-independent methods can be found in Nesseris:2013PhRvD . PCA is a multivariate analysis and is usually employed to predict the form of cosmological quantities in a model-independent, non-parametric manner Huterer2003prl ; clarkson2010prl ; Huterer2005 ; Zheng:2017ulu ; Ishida2011is . In Qin:2015eda and Liu:2015yha , different variants of PCA techniques have been adopted. Reference Qin:2015eda uses an error model and then creates different set of simulated Hubble data to construct the covariance matrix while Liu:2015yha use the weighted least square method and combines it with PCA.

From the viewpoint of the application of PCA, there exist two distinct methodologies. These methods differ mainly in the way the covariance matrix is calculated, which is the first step of any PCA technique. One way to implement PCA is through the computation of Fisher matrix, Nesseris:2013PhRvD ; Huterer2003prl ; clarkson2010prl ; Huterer2005 ; Zheng:2017ulu ; Ishida2011is ; Crittenden:2005wj ; Miranda:2018 ; Hart:2019 ; Hojjati:2012prd ; 2013JCAP…02..049N ; Hart:2021kad . One can bin the redshift range and assume a constant value for the quantity to be reconstructed in that redshift bin. These constant values are the initial parameters of the PCA. Fisher Matrix quantifies the correlation and uncertainties of these parameters. Therefore, by deriving these constants in different bins using PCA, we can reproduce our targeted quantity in terms of redshift Ishida2011is . Alternatively, a polynomial expression for the dynamical quantity can be assumed. In this case, the coefficients of the polynomial are the initial parameters for PCA, and the analysis gives the final values of these coefficients, which eventually gives the dynamical quantity clarkson2010prl .

PCA is independent of any prior biases and it is also helpful in comparing the quality of different data-sets Crittenden:2005wj ; yu2013prd . It is an application of linear algebra, which makes the linearly correlated data points uncorrelated. The correlated data points of the data-set used in PCA are transformed by rotating the axes, where the angle of rotation of these axes is such that linear-correlations between data-points are the smallest compared to any other orientation. The new axes are the principal component(PC)s of the data points, and these PCs are orthogonal to each other. In terms of information in the data-set, PCA creates a hierarchy of priority between these PCs. The first PC contains information of the signal the most and hence has the smallest dispersion of data-points about it. The second PC contains less information than the first PC and therefore provides a higher dispersion of the data-points as compared to the first PC. Higher-order PCs have the least priority as these correspond to noise and we can drop them. The reduction of dimensions is a distinctive feature of PCA. Therefore the final reconstructed curve in the lower dimension corresponds predominantly to the signal of the data-set. The PCA method also differs from the regression algorithms, which can not distinguish between signal and the noise. PCA can omit the features coming from the noise part and can pick the actual trend of the data-points licia2010lnip ; Sivia2006 ; Steinhardt2018 . In our case, the data-points on which we apply PCA are created in the process of reconstruction, see sect(2).

We assume polynomial expressions for the observational quantities, namely the Hubble parameter $H(z)$, and the distance modulus ($\mu(z)$). The polynomial form is then modified and predicted solely by PCA reconstruction algorithm. The only assumption we make is that it is possible to expand the reconstructed quantity in a polynomial. The dark energy equation of state is constructed in the subsequent steps after the reconstruction of the Hubble parameter and distance modulus.

The two-step process of reconstruction is termed as the derived approach. Further, in a direct approach, we apply the PCA technique to the polynomial expression of the equation of state parameter $w(z)$. For the application of our methodology, we require a tabulated dataset of the observational quantity, which can schematically be represented as, (independent variable)–(dependent variable)–(error bar of dependent variable). In the direct approach, described below, we do not have the tabulated data-set of $w(z)$ and we derive the functional form of $w(z)$ from the tabulated data-set of $z-H(z)$ and $z-\mu(z)$. In the case of the derived approach, PCA first reconstructs the functional form of $H(z)$ and $\mu(z)$ from the Hubble parameter and Supernovae data-set respectively without any cosmological model. In the absence of cosmological model, $H(z)$ and $\mu(z)$ have no any parameter dependencies. However to construct the final form of equation of state of dark energy, $w(z)$ we need $\Omega_{m}$ as input. Application of PCA on polynomial expressions to reconstruct cosmological quantities have been employed earlier in Ishida2011is ; Qin:2015eda ; Liu:2015yha .

For a further check on the robustness of our results; we compute the correlations of the coefficients, present in the polynomial expressions of the reconstructed quantities. We apply PCA in the coefficient space, which is created by the polynomial. The PCA rearranges the correlation of these coefficients. Calculation of correlation-coefficients help to identify the linear and non-linear correlations of the components in the reconstructed quantities. We show that using the correlation coefficients calculation, we can restrict the allowed terms in the polynomial.

This paper is structured as follows. In sect(2) we describe the reconstruction algorithm and the use of correlation tests in those reconstructions. In sect(3), we describe the two approaches of reconstruction we follow. We present our results in sect(4). In sect(5), we conclude by summarizing the results.

We begin with an initial basis, $g_{i}=f(x)^{i-1}$, where $i=1,2,....,N$ through which we can express the quantity to be reconstructed as,

The initial basis can be written in matrix form as, $\mathbf{G}=(f_{1}(x),f_{2}(x),...,f_{N}(x))$. Coefficients ${b_{i}}$ create a coefficient space of dimension $N$. Each point in the coefficient space gives a realization of $\xi(x)$, including a constant $\xi(x)$. PCA modifies the coefficient space and chooses a single realization. One can consider different kinds of functions as well as different combinations of polynomials as initial bases. Choosing one initial basis function over the other is done by correlation coefficient calculation described below.

The Pearson correlation coefficient for two parameters $A$ and $B$ is given by,

where $\rho\in[-1,1]$. For linearly uncorrelated variables, the correlation coefficient, $\rho=0$. An exact correlation is identified by $\rho=-1$ or $\rho=+1$. The Spearman rank coefficient is, in turn, the Pearson correlation coefficient of the ranks of the parameters; rank being the value assigned to a set of objects and it determines the relation of every object in the set with the rest of them. We mark the highest numeric value of a variable $A$ as ranked 1, the second-highest numeric value of the variable as ranked 2 and so on. A similar ranking is done for the ranks of the parameter $B$.

To obtain the coefficients for the correlation analysis, we divide the parameter space into $n$ patches. We therefore have $n$ values associated with one coefficient of the polynomial; this is the number of columns of the coefficients matrix $\mathbf{Y}$, (eq(5)). After ranking all the values of $A$ and $B$, we obtain the table for the ranks of $A$ and $B$. We then proceed to compute the Spearman correlation coefficient($r$), which is the Pearson Correlation coefficient of rank of $A$ and $B$. Like in the case of the Pearson Correlation coefficients, $r\in[-1,1]$. Computing Kendall correlation coefficient($\tau$) gives a prescription to calculate the total number of concordant and dis-concordant pairs from the values of the variables $A$ and $B$ M_G_Kendall_1938 .

If we pick two pairs of points from the table of $A$ and $B$, say $(a_{i},b_{i})$ and $(a_{j},b_{j})$, for $i\neq j$ if $a_{i}>a_{j}$ when $b_{i}>b_{j}$ or if $a_{i}<a_{j}$ when $b_{i}<b_{j}$; then that pair of points are said to be in concordance with each other. On the other hand, for $i\neq j$, $a_{i}>a_{j}$ when $b_{i}<b_{j}$ or if $a_{i}<a_{j}$ when $b_{i}>b_{j}$, then these two pairs are called to be in dis-concordance with each other. Every concordant pair is scored $+1$ and every dis-concordant pair is scored $-1$. The Kendall correlation coefficients are defined as,

Again, if $N_{cp}$ is the number of concordance pair and $N_{dp}$ is the number of dis-concordance pair

Hence the expression of $\tau$ is,

where, $\tau\in[-1,1]$.

We perform the correlation coefficients calculation twice. The first time is to select the number of terms in the initial polynomial $N$(eq(1)), and second time to select the number of terms in the final polynomial, (eq(8)). We select that value of $N$ for which the Pearson Correlation Coefficient is higher than the Spearman and Kendall Correlation coefficients. We use the R-package for statistical computing to calculate the correlation coefficients R:manual .

Values of linear and non-linear correlations depend on the quantity we want to reconstruct. They are also sensitive towards the data-set that we use in reconstruction. For instance, reconstruction of a fast varying function, which has non-zero higher order derivatives will introduce more non-linear contributions to the correlation of $b_{i}$ than linear contributions, and we need a greater set of initial basis, which means a higher value of $N$. We take different values of $N$ and check the linear and non-linear correlation coefficients to fix the value of $N$. Though a large value would help, we can not, however, fix $N$ to any arbitrarily significant number as it makes the analysis computationally expensive.

To compute the covariance matrix, we define the coefficient matrix ($\mathbf{Y}$) by selecting different patches from the coefficient space,

where $n$ is the number of patches that we have taken into account and $N$ is the total number of initial basis defined in eq(1); therefore, Y is a matrix of dimension $N\times n$ and $b_{n}^{(N)}$ being the value of Nth coefficient in $n$th patch.

In the present analysis, we have taken $n$ to be the order of $10^{3}$. We estimate the best-fit values of the coefficients at each patch by $\chi^{2}$ minimization, where $\chi^{2}$ is defined as

$k$ is the total number of points in the data-sets. If the observational data-set have significant non-diagonal elements in the data covariance matrix($C_{data}$), we have to incorporate $C_{data}$ in eq(6), as $\chi^{2}=\Delta^{T}\mathcal{C}_{data}^{-1}\Delta$, where $\Delta=\xi(x)_{data}-\xi(\{b_{i}\},x)$ and $x$, $\xi(x)_{data}$ vary for each data-points. Calculation of $\chi^{2}$ in all the patches gives us the variation of the $N$ coefficients and finally gives $n$ number of points in the coefficient space, over which we apply PCA. We calculate the covariance matrix and correlations of the coefficients for these $n$ points. In this analysis, each patch contains the origin of the multi-dimensional coefficient space.

The covariance matrix of the coefficients, $\mathbf{C}$ is written as,

Eigenvector matrix, $\mathcal{E}$ of this covariance matrix will rotate the initial basis of the coefficient space to a position where the patch-points will be uncorrelated. We organize the eigenvectors in the eigenvector matrix $\mathcal{E}$ in the increasing order of eigenvalues. Eigenvalues of the Covariance matrix quantifies the error associated with each principal component Huterer2003prl ; Ishida2011is ; 2019Gortler ; Louis1953 .

If $\mathbf{U}=(u_{1}(x),u_{2}(x),...,u_{N}(x))$ the final basis is given by,

The final reconstructed form of $\xi(x)$ is,

where $M\leq N$ and the $\beta_{i}$s are the uncorrelated coefficients associated with the final basis. The coefficients $\beta_{i}$ are re-calculated by using $\chi^{2}$ minimization, for the same $n$ patches considered earlier to create initial coefficient matrix (5). The value of $M$ can be determined by correlation coefficient calculation discussed above. As PCA only breaks the linear correlation, we select $M$ for which PCA is able to break Pearson correlation coefficient to the largest extent. Eigenvalues of the covariance matrix ($e_{1},e_{2},e_{3},...,e_{M}$) indicate how well PCA can pick the best patch-point in the coefficient space. We look for the lowest number of final basis which can break the linear correlation to the greatest extent. We can start with smaller value of $N$ which eventually influence the value of $M$, but it will pose the risk of losing essential features from PCA data-set.

The Hubble parameter ($H(z)$) for a spatially flat Universe, composed of dark energy and non-relativistic matter is given by,

Here we have assumed that the contributions to the energy density is only due to the non-relativistic dark matter and dark energy. The density parameters for non-relativistic matter and dark energy are given by $\Omega_{m}$ and $\Omega_{x}$. The quantity $H_{0}$ denotes the present-day value of the Hubble parameter, namely the Hubble constant and $w(z)$ is the dark energy equation of state parameter ($w(z)$).

We assume no interaction between matter and dark energy in the present analysis. In the following subsections, we discuss the derived and direct approach to reconstruct $w(z)$ using PCA.

In this paper, we reconstruct late-time cosmology using the Principal Component Analysis. There are very few prior assumptions about nature and distribution of different components contributing to the energy of the Universe. Observational Hubble parameter and distance modulus measurements of type Ia supernovae are the observable quantities that are taken into account in the present analysis. We proceed in two different ways to do the reconstruction. The first one is a derived approach where the observable quantities are reconstructed from the data using PCA, and then $w(z)$ is obtained from the reconstructed quantities using Friedman equation. The other approach is a direct one. In this case, $w(z)$ is reconstructed directly from the observational data using PCA without any intermediate reconstruction. Based on the efficiency of the method to break the correlation among the coefficients we can select one reconstruction curve over the other. We achieve a better reconstruction in the derived approach as compared to the direct approach, even though the direct approach has lesser uncertainties than the derived approach in predicting the patch of the best coefficient point from the $N$ dimensional coefficient space due to the lower number of initial and final bases. In the derived approach, though, we need no cosmological theory or model to determine the functional form of $H(z)$ and $\mu(z)$; however for the calculation of $w(z)$, we need to input a cosmological model. On the other hand, the direct approach starts with using a polynomial form of $w(z)$ in Friedman equation. Therefore, in the case of direct approach, the dependencies of $w(z)$ on the cosmological parameters are greater.

We have adopted the simulated as well as observed data-sets for our analysis. Simulated data-sets are used to check the efficiency. For the reconstruction of $w(z)$ the analysis produces consistent result only for the Hubble parameter data. Though the reconstruction of $\mu(z)$ through the derived approach is consistent and reconstruction of $\mu(z)$ is within the error bars of the data-set, the result for the reconstructed $w(z)$ deviates drastically from the physically acceptable range. The increase in the order of differentiation to connect $w(z)$ with $\mu(z)$ is a possible reason for this inconsistency. Here we are focusing on the $w(z)$ reconstruction only from $H(z)-z$ data-set. The reconstructed $w(z)$ by the variable $(1-a)$, obtained in the derived approach for $H(z)-z$ data-set shows a phantom like nature, that is $w(z)<-1$ at present and a non-phantom nature in the past for most of the values of $\Omega_{m}$. In the case of direct approach, $w(z)$ curves show oscillations in the phantom and non-phantom regime. The calculation of the correlation coefficients clearly shows a preference for the derived approach. PCA lacks the efficiency of breaking the correlation in the initial basis in case of the direct approach. This probably causes the inconsistency between the results obtained in the derived and direct approaches.

The other important factor is the variable of reconstruction. In the present analysis, we have adopted three different reconstruction variables, namely $(1-a)$, $a$ and $z$, in both direct and derived approaches. Values of correlation coefficients after PCA select the reconstruction variable $(1-a)$ over the other two. We should emphasize the result obtained for variable $(1-a)$ by derived approach. Since we have a finite number of terms in the initial polynomial, due to the constraints set by computational power, to some extent the results depends on the assumption of polynomial expression of the observables. One of our future plans is to develop an algorithm which can inclusively select the most suitable initial basis form for a fix computational power and a given observational data-set. The reconstructed curves, obtained for $(1-a)$, show that $w(z)$ shows a phantom nature at present epoch and it was in non-phantom nature in the past. PCA reconstruction indicates an evolution of the dark energy equation of state parameter.

The correlation coefficients of the first three coefficients of the polynomial expression eq(1) and eq(8) of the reconstructed quantity is shown in table 3. PCA breaks the linear correlation of the coefficients, which is evident from table[3]. Presence of the non-linear correlation in the initial coefficients complicates the process. As the first three terms of the ultimate expression of the reconstructed quantity contains the dominant trend of the data-points, we only mention the correlation coefficients only for the first three parameters. The reconstruction which can break Pearson Correlation to a greater extent as well as have lesser Spearman and Kendall Correlation coefficients is selected. We can see from table 3 that the reconstruction by $(1-a)$, breaks the correlation to a greater extent than in the case of $a$ and $z$. Variation of correlation coefficients before and after the application of PCA is similar for both the simulated as well as the real data-set. Difference of Pearson correlation coefficients for simulated and real data-set in case of $(1-a)$ and $a$ is of the order of $10^{-6}$ or less. For the variable $z$, the difference of Pearson correlation coefficients for real and simulated data is of the order of $10^{-4}$ or less.

The authors acknowledge the use of High-Performance-Computing facility at IISER Mohali. AM acknowledges the financial support from the Science and Engineering Research Board (SERB), Department of Science and Technology, Government of India as a National Post-Doctoral Fellow (NPDF, File no. PDF/2018/001859). The authors would like to thank J S Bagla for useful suggestions and discussions.

The observational data-sets used in the analysis are publicly available and duly referred in the text. The simulated data-sets can be provided on request with appropriate justification.