We know that if we have a function with several variables $x_1,\dots, x_n$, the second degree Taylor expansion around $a=(a_1,\dots,a_n)$ gives
$f(x_1, ... , x_n) \approx f(a)+Df(a)(x-a)+\frac{1}{2}(x-a)^THf(a)(x-a), $
where $x=(x_1,\dots,x_n)$, $Df(x)$ represents the $1\times n$ gradient matrix and $Hf(x)$ represents the $n\times n$ hessian matrix.
Here, each variable is real-valued, i.e., $x_i \in \mathbb{R}, i=1,\dots,n.$
But suppose now that each variable is a vector variable, i.e., $\boldsymbol{x}_i \in \mathbb{R}^p, i=1,\dots,n.$ Does the Taylor formula be the same for $f(\boldsymbol{x}_1, ... , \boldsymbol{x}_n)$ around $a=(\boldsymbol{a}_1,\dots,\boldsymbol{a}_n)$?