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 . We demonstrate the applicability of our results by applying them 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.