Gaussian Process: vector valued response Gaussian Processes (GPs) define a prior over functions that can be updated to a posterior once we have observed data. I've been working with scalar-valued GPs, i.e. functions $f: \mathbb{R}^{d} \rightarrow \mathbb{R}$ such that
\begin{equation}
f(x_1,...,x_d) \sim GP(m(x), K(x,x'))
\end{equation}
However, I was wondering if Gaussian Processes can be used for vector-valued response, i.e. when $f: \mathbb{R}^{d} \rightarrow \mathbb{R}^{n}$? This could be useful if you are trying to jointly estimate $n$ parameters.
 A: Yes, Gaussian processes can be used to approximate vector valued functions. 
This is done in spatial statistics where it's called cokriging and in machine learning in the field of multi-task learning. There's a relatively recent survey paper by Alvarez et al. which includes Gaussian processes this as part of a wider survey of kernel methods for vector valued functions. There are about ~120 reference spread across the ML and stats literature so the topic has attraced a fair amount of academic interest. 
For comparison to what you've discussed above, I think the definition given by Alvarez et al. was instructive. 
Assume you're interested in learning a function $f$ where $f: \chi = \mathbb{R}^p \rightarrow \mathbb{R}^D$, and you have training data $S = (\mathbf{X}, \mathbf{Y}) = \{(x_1, y_1), . . . ,(x_N , y_N )\}$. Then:

The vector-valued function $f$ is assumed to follow a Gaussian process
  $$ f ∼ GP(\mathbf{m}, \mathbf{K}),  $$ where $\mathbf{m} \in
 \mathbb{R}^D$ is a vector which components are the mean functions
  ${m_d(x)}^D_{d=1}$ of each output and $\mathbf{K}$ is a positive
  matrix valued function [that is, $\mathbf{K} : \chi \times \chi
 \rightarrow R^{D \times D}$, such that for any $\mathbf{x}, \mathbf{x
 ′}$ we have $\mathbf{K}(\mathbf{x}, \mathbf{x}′)$ is a positive
  semi-definite matrix]. The entries  $(\mathbf{K}(\mathbf{x},
 \mathbf{x}′))_{d,d′}$ in the matrix $\mathbf{K}(\mathbf{x},
 \mathbf{x}′)$ correspond to the covariances between the outputs
  $f_d(\mathbf{x})$ and $f_{d′} (\mathbf{x} ′)$ and express the degree
  of correlation or similarity between them. For a set of inputs
  $\mathbf{X}$, the prior distribution over the vector $f(\mathbf{X})$
  is given by $$ f(\mathbf{X}) ∼ N (m(\mathbf{X}), \mathbf{K}(\mathbf{X},
 \mathbf{X})), $$ where $m(\mathbf{X})$ is a vector that concatenates
  the mean vectors associated to the outputs and the covariance matrix
  $\mathbf{K}(\mathbf{X}, \mathbf{X})$ is the block partitioned matrix
  where $m(\mathbf{X})$ is a vector that concatenates the mean vectors
  associated to the outputs and the covariance matrix
  $\mathbf{K}(\mathbf{X}, \mathbf{X})$ is the block partitioned matrix
  [given by the $ND \times ND$ entries $(\mathbf{K}(\mathbf{x}_i,
\mathbf{x}_j ))_{d,d′}$  for $i,j\in \{1 \ldots N \}$ and $d, d' \in \{1 .. D \}$]

Alvarez, Mauricio A., Lorenzo Rosasco, and Neil D. Lawrence. "Kernels for vector-valued functions: A review." Foundations and Trends® in Machine Learning 4.3 (2012): 195-266.
I do recommend the paper, I thought it was very interesting  and I hope
that helps!
