In this paper we develop a framework for multivariate functional approximation by a suitable Gaussian process via an exchangeable pairs coupling that satisfies a suitable approximate linear regression property, thereby building on work by Barbour (1990) and Kasprzak (2020). We demonstrate the applicability of our results by applying it to joint subgraph counts in an Erdős-Renyi random graph model on the one hand and to vectors of weighted, degenerate $U$-processes on the other hand. As a concrete instance of the latter class of examples, we provide a bound for the functional approximation of a vector of success runs of different lengths by a suitable Gaussian process which, even in the situation of just a single run, would be outside the scope of the existing theory.