  In the current work, we have implemented an extension of the standard
Gaussian Process formalism, namely the Multi-Task Gaussian Process with the
ability to perform a joint learning of several cosmological data
simultaneously. We have utilised the "low-redshift" expansion rate data from
Supernovae Type-Ia (SN), Baryon Acoustic Oscillations (BAO) and Cosmic
Chronometers (CC) data in a joint analysis. We have tested several possible
models of covariance functions and find very consistent estimates for
cosmologically relevant parameters. In the current formalism, we also find
provisions for heuristic arguments which allow us to select the best-suited
kernel for the reconstruction of expansion rate data. We also utilised our
method to account for systematics in CC data and find an estimate of $H_0 =
68.52^{+0.94 + 2.51 (sys)}_{-0.94} $ $\textrm{km/s Mpc}^{-1}$ and a
corresponding $r_d = 145.61^{+2.82}_{ - 2.82 - 4.3 (sys)} $ Mpc as our primary
result. Subsequently, we find constraints on the present deceleration parameter
$q_0 = -0.52 \pm 0.06$ and the transition redshift $z_T =
0.64^{+0.12}_{-0.09}$. All the estimated cosmological parameters are found to
be in good agreement with the standard $\Lambda$CDM scenario. Including the
local model-independent $H_0$ estimate to the analysis we find $H_0 = 71.40^{ +
0.30 + 1.65 (sys)}_{- 0.30 } $ $\textrm{km/s Mpc}^{-1}$ and the corresponding
$r_d = 141.29^{ + 1.31 }_{-1.31-2.63 (sys)}$ Mpc. Also, the constraints on $r_d
H_0$ remain consistent throughout our analysis and also with the
model-dependent CMB estimate. Using the $\mathcal{O}m(z)$ diagnostic, we find
that the concordance model is very consistent within the redshift range $z
\lesssim 2$ and mildly discrepant for $z \gtrsim 2$.
