Let $f:\mathbb{R}^m \to \mathbb{R}^n$. Is it possible to do a Taylor expansion of $f$ around $\theta\in\mathbb{R}^m$? I am hoping for something like

$$f\left(\theta\right) = f\left(\theta_0\right) + A \left(\theta - \theta_0\right) + \left(\theta - \theta'\right)^T \text{something} \left(\theta - \theta'\right)$$

where $A$ is a $n\times m$ matrix, and its rows are gradient of $f_i$ ($i$-the entry of vector $f$) with respect to vector $\theta$. But what should "something" be?

For example, let us consider the simple case where $f :\mathbb{R}^m \to \mathbb{R}^2$. Denote $f = \left(f_1,f_2\right)$. Then

$$f_1\left(\theta_n\right) = f_1(\theta_0) + \nabla f_1\left(\theta_0\right)(\theta_n - \theta_0) + (\theta_n - \theta_0)^T H_1(\theta') (\theta_n - \theta_0) $$

$$f_2\left(\theta_n\right) = f_2(\theta_0) + \nabla f_2\left(\theta_0\right)(\theta_n - \theta_0) + (\theta_n - \theta_0)^T H_2(\theta') (\theta_n - \theta_0) $$

If I try to put this in matrix form, I will get something like

$$f\left(\theta_0\right) = f(\theta_0) + \begin{pmatrix} \nabla f_1(\theta_0)\\ \nabla f_2(\theta_0)\end{pmatrix} (\theta_n - \theta_0) + ???$$