Unfortunately, Lagrange-type results do not hold for vector-valued functions. In other words, it is not true that there exists an $\tilde{\x}$ such that
\[\f(\x) = \f(\x_0) + \nabla \f(\tilde{\x}) \tp (\x - \x_0);\]such a point exists for each element of $\f$ separately, but these points will not be the same. The following theorem states this explicitly.
Theorem: Suppose $\f:\real^d \to \real^k$ is differentiable on $N_r(\x_0)$. Then for any $\x \in N_r(\x_0)$,
\[\f(\x) = \f(\x_0) + \D(\x, \x_0) (\x - \x_0),\]where row $j$ of $\D(\x, \x_0)$ is $\nabla f_j(\tilde{\x}_j) \tp$, with each $\tilde{\x}_j$ on the line segment between $\x$ and $\x_0$.
Proof: The above equation is simply repeating the Lagrange Taylor expansion for each component, then stacking the results into a linear system of equations.
This result can be simplified if the second partial derivatives of $\f$ are bounded on $N_r(\x_0)$: $\abs{\partial^2 f_j / \partial x_m \partial x_l} \le M$ for all $j, m, l$. In this case,
\[\D(\x, \x_0) = \nabla \f(\x_0) \tp + \R, \qquad \norm{\R}_{\max} \le \sqrt{d} M \norm{\x - \x_0},\]where $\norm{\R}_{\max}$ is the largest element of $\R$ in absolute value.
Proof: Element $(j,m)$ of $\R$ is $\partial_m f_j(\tilde{\x}_j) - \partial_m f_j(\x_0)$, where $\partial_m f_j = \partial f_j / \partial x_m$ is a scalar-valued function. Its gradient has $d$ elements, each bounded by $M$, so its norm is at most $\sqrt{d} M$. Thus, for some $\x^\ast$ between $\tilde{\x}_j$ and $\x_0$,
\[\begin{alignat*}{2} \abs{\partial_m f_j(\tilde{\x}_j) - \partial_m f_j(\x_0)} &= \abs{\nabla (\partial_m f_j)(\x^\ast) \tp (\tilde{\x}_j - \x_0)} &\hspace{4em}& \href{taylor-expansion.html}{\text{Mean value theorem}} \\ &\le \norm{\nabla (\partial_m f_j)(\x^\ast)} \norm{\tilde{\x}_j - \x_0} && \href{cauchy-schwarz.html}{\text{Cauchy-Schwarz}} \\ &\le \sqrt{d} M \norm{\x - \x_0} && \tilde{\x}_j \text{ lies between } \x \text{ and } \x_0 \end{alignat*}\]In statistics, this theorem is usually applied with respect to the score.