Submodular functions have found a wealth of new applications in data science and machine learning models in recent years. This has been coupled with many algorithmic advances in the area of submodular optimization: (SO) $\min/\max~f(S): S \in \mathcal{F}$, where $\mathcal{F}$ is a given family of feasible sets over a ground set $V$ and $f:2^V \rightarrow \mathbb{R}$ is submodular. In this work we focus on a more general class of \emph{multivariate submodular optimization} (MVSO) problems: $\min/\max~f (S_1,S_2,\ldots,S_k): S_1 \uplus S_2 \uplus \cdots \uplus S_k \in \mathcal{F}$. Here we use $\uplus$ to denote disjoint union and hence this model is attractive where resources are being allocated across $k$ agents, who share a `joint' multivariate nonnegative objective $f(S_1,S_2,\ldots,S_k)$ that captures some type of submodularity (i.e. diminishing returns) property. We provide some explicit examples and potential applications for this new framework. For maximization, we show that practical algorithms such as accelerated greedy variants and distributed algorithms achieve good approximation guarantees for very general families (such as matroids and $p$-systems). For arbitrary families, we show that monotone (resp. nonmonotone) MVSO admits an $α(1-1/e)$ (resp. $α\cdot 0.385$) approximation whenever monotone (resp. nonmonotone) SO admits an $α$-approximation over the multilinear formulation. This substantially expands the family of tractable models for submodular maximization. For minimization, we show that if SO admits a $β$-approximation over \emph{modular} functions, then MVSO admits a $\frac{β\cdot n}{1+(n-1)(1-c)}$-approximation where $c\in [0,1]$ denotes the curvature of $f$, and this is essentially tight. Finally, we prove that MVSO has an $αk$-approximation whenever SO admits an $α$-approximation over the convex formulation.