Bivariate partial information decompositions (PIDs) characterize how the information in a "message" random variable is decomposed between two "constituent" random variables in terms of unique, redundant and synergistic information components. These components are a function of the joint distribution of the three variables, and are typically defined using an optimization over the space of all possible joint distributions. This makes it computationally challenging to compute PIDs in practice and restricts their use to low-dimensional random vectors. To ease this burden, we consider the case of jointly Gaussian random vectors in this paper. This case was previously examined by Barrett (2015), who showed that certain operationally well-motivated PIDs reduce to a closed form expression for scalar messages. Here, we show that Barrett's result does not extend to vector messages in general, and characterize the set of multivariate Gaussian distributions that reduce to closed-form. Then, for all other multivariate Gaussian distributions, we propose a convex optimization framework for approximately computing a specific PID definition based on the statistical concept of deficiency. Using simplifying assumptions specific to the Gaussian case, we provide an efficient algorithm to approximately compute the bivariate PID for multivariate Gaussian variables with tens or even hundreds of dimensions. We also theoretically and empirically justify the goodness of this approximation.