Gaussian process regression (GPR) is a Bayesian non-parametric technology that has gained extensive application in data-based modelling of various systems, including those of interest to chemometrics. However, most GPR implementations model only a single response variable, due to the difficulty in the formulation of covariance function for correlated multiple response variables, which describes not only the correlation between data points, but also the correlation between responses. In the paper we propose a direct formulation of the covariance function for multi-response GPR, based on the idea that its covariance function is assumed to be the “nominal” uni-output covariance multiplied by the covariances between different outputs. The effectiveness of the proposed multi-response GPR method is illustrated through numerical examples and response surface modelling of a catalytic reaction process.