Multivariate Taylor Expansion I am in confidence with Taylor expansion of function $f\colon R \to R$, but I when my professor started to use higher order derivatives and multivariate Taylor expansion of $f\colon R^n \to R$ and $f\colon R^n \to R^m$ I felt lost.
Can somean explain to me from scratch multivariate Taylor?
In particular I don't understand the notation
$$
f(x+h) = \sum_{k=0}^p \frac{1}{k!} f^{(k)}(x)[h,...,h] + O(h^{p+1})
$$
Why we need the k-linear form $ \frac{1}{k!} f^{(k)}(x)[h,...,h]$? This k-linear form is the derivative or the derivative is only $f^{(k)}(x)$?
I'm quite lost. Thank you.
 A: It is enough to understand the case $f:\ {\mathbb R}^n\to {\mathbb R}$, and we expand $f$ at $x=0$. So we are interested in polynomials of low degree $d$ in the variables $x_1$, $\ldots$, $x_n$ that approximate $f$ in the neighborhood of $0\in{\mathbb R}^n$.
The best polynomial of degree $\leq0$ for this purpose is obviously the constant function $$j^0:\quad x\mapsto j^0(x)= f(0)\ .$$ In fact we have
$$f(x)=j^0(x)+ r(x),\quad \lim_{x\to0} r(x)=0\ .$$
The best approximating polynomial of degree $\leq1$ is the  function $$j^1:\ x\mapsto j^1(x)=f(0)+\nabla f(0)\cdot x\ ,$$
where $\nabla f(0):=\left({\partial f\over\partial x_1},\ldots,{\partial f\over\partial x_n}\right)_0$.
In fact we have
$$f(x)=j^1(x)+ r(x),\quad \lim_{x\to0} {|r(x)|\over |x|}=0\ .$$
Advancing to $d=2$ we have to be aware that the dimension of the space of homogeneous polynomials of degree $2$ is ${n(n+1)\over2}$, because we have to envisage all terms of the form $x_k^2$ and also mixed terms $x_ix_j$ with $i\ne j$. It follows that $j^2(x)$ will have the already known terms of degree $0$ and $1$ and "quadratically in $n$ many" terms of degree $2$. The latter are best arranged in a matricial way, so that one has a sum $\sum_{1\leq i\leq n,\ 1\leq j\leq n}a_{ij} x_ix_j$. The coefficients $a_{ij}$ will be related to the second partials of $f$ at $0$ in a particular way established in the proof of the general theorem. Going into the details one obtains
$$j^2(x)=f(0)+\nabla f(0)\cdot x+{1\over2}\sum_{1\leq i\leq n,\ 1\leq j\leq n}\left({\partial^2 f\over\partial x_i\partial x_j}\right)_0 x_ix_j\ ,$$
and one can prove that
$$f(x)=j^2(x)+ r(x),\quad \lim_{x\to0} {|r(x)|\over |x|^2}=0\ .$$
The formula you displayed is a systematic and condensed way to write this all up. E.g., in degree $3$ we shall obtain a triple sum of $n^3$  monomials $x_i x_j x_k$, each of them multiplied with the appropriate third partial derivative of $f$ at the origin.
A: One can think about Taylor's theorem in calculus as applying in the following cases:


*

*Scalar-valued functions of a scalar variable, i.e. $f : \mathbb{R} \rightarrow \mathbb{R}$

*Vector-valued functions of a scalar variable, i.e. $f : \mathbb{R} \rightarrow \mathbb{R}^n$

*Scalar-valued functions of a vector variable, i.e. $f : \mathbb{R}^n \rightarrow \mathbb{R}$

*Vector-valued functions of a vector variable, i.e. $f : \mathbb{R}^n \rightarrow \mathbb{R}^m$


All of these can be derived & proven based on nothing more than integration by parts (the last one needs to be developed in a banach space & the third one is more commonly reduced to the first one which is just a shorthand for re-proving it via integration by parts) if you set things up correctly (as is done in Lang's Undergraduate, Real & Functional Analysis books) & so your main obstacle here is formalism - this is no small obstacle as we'll see below. 
Now I'm not sure if your expression for Taylor's formula is map 3 or map 4, one would think it is map 3 since you used the word "linear form" which is standard parlance for maps from vector spaces into a field but you did ask about maps of the form $f : \mathbb{R}^n \rightarrow \mathbb{R}^m$ - are you sure you are differentiating these kinds of maps because they add a whole world of complexity compared to the first 3? 
If you're asking about maps of the form in 3 then some intuition is given in this video & some examples & a proof are given in this video. After these you should have enough of a grasp of what's going on & if you focus on developing the formalism properly you should be able to prove it yourself in more general spaces.
If you're actually asking about map 4 then you may be used to the definition of the derivative of the last map as $f : \mathbb{R}^n \rightarrow \mathbb{R}^m$ as something like $\mathcal{f'} : \mathbb{R}^n \rightarrow \mathcal{L}(\mathbb{R}^n,\mathbb{R}^m)$ satisfying all the conditions a differentiable map does, well it's second derivative is defined similarly using a map of the form $f'': \mathbb{R}^n \rightarrow \mathcal{L}(\mathbb{R}^n,\mathcal{L}(\mathbb{R}^n,\mathbb{R}^m))$ & so on, however you see no such entity as a "k-linear form" in any of this & that's because there is a theorem which allows one to think of maps like the second derivative above in terms of multilinear maps & so one can re-cast the theory using multilinear maps which eases the development & allows for nice proofs etc... but without this being explained it might appear odd to randomly start pulling out multilinear maps.
In any case the derivative is a linear map by definition & so that is why you're coming across the word linear, but since you didn't put subscripts on the $[h,...,h]$ I'm not sure how deep I can go, because I see two possibilities here so if the above isn't enough of an explanation as it stands just let me know.
A: For 2 variables I believe the following expression to be correct. If it is incorrect then I should like to know!
$$f(x,y) = \sum_{m=0}^{\infty} \sum_{n=0}^{\infty} \beta_{m,n} (x-x_0)^m (y-y_0)^n$$
$$\beta_{m,n} = \frac{1}{n! m!} \left[ \frac{\partial^{m+n}}{\partial x^m \partial y^n} \ f(x,y) \right]_{(x=x_0,y=y_0)}$$
Where the square bracket notation with coordinate subscript means evaluate at.
I previously derived the following formula, however I now suspect that it is not correct. Perhaps someone better than myself at proof by induction could verify this?
$$f(x,y) = \sum_{k=0}^{\infty} \frac{1}{k!} \sum_{j=0}^{k} \alpha_{k,j} (x-x_0)^{j} (y-y_0)^{k-j}$$
$$\alpha_{k,j} = \,^{k}C_{j} \left[ \frac{\partial^{k}}{\partial x^j \partial y^{k-j}} \ f(x,y) \right]_{(x=x_0,y=y_0)}$$
In order to prove the first form, I did the following:
$$\left[ f(x,y) \right]_{(x0,y0)} = \sum_{m=0}^{\infty} \sum_{n=0}^{\infty} \beta_{m,n} \cdot 0^m \cdot 0^n = \beta_{0,0}$$
Differentiating partially with respect to x we find
$$\left[ \frac{\partial f(x,y)}{\partial x} \right]_{(x0,y0)} = \sum_{m=0}^{\infty} \sum_{n=0}^{\infty} \left( \beta_{m,n} \cdot m \cdot 0^{m-1} \cdot 0^n + \left(\frac{\partial \beta_{m,n}}{\partial x} \right) \cdot0^m\cdot0^n\right) = \beta_{1,0} $$ where $\tfrac{\partial \beta_{m,n}}{\partial x} = 0$.
Continuing this process for all combinations of partial derivatives with respect to both x and y, we find that in general $$\beta_{m,n}=\frac{1}{m!n!}\cdot \left[ \frac{\partial f(x,y)}{\partial x} \right]_{(x=x0,y=y0)}$$
Unless I have made a mistake that is. (In which case please add a comment so I can correct my answer.)
A: For 3 variables:
$$f(x,y,z)=f(x_0,y_0,z_0)$$ 
$$+\frac{\partial f_0}{\partial x}(x-x_0)+\frac{\partial f_0}{\partial y}(y-y_0)+\frac{\partial f_0}{\partial z}(z-z_0)\quad \Rightarrow Order 1$$
$$+\frac{1}{2} \bigg(\frac{\partial^2 f_0}{\partial x^2}(x-x_0)^2+\frac{\partial^2 f_0}{\partial y^2}(y-y_0)^2+\frac{\partial^2 f_0}{\partial z^2}(z-z_0)^2+2\frac{\partial^2 f_0}{\partial x\partial y}(x-x_0)(y-y_0) $$
$$+2\frac{\partial^2 f_0}{\partial x\partial z}(x-x_0)(z-z_0)+2\frac{\partial^2 f_0}{\partial z\partial y}(z-z_0)(y-y_0)\bigg)\quad \Rightarrow Order 2$$
And it goes like this to higher orders
A: Let's suppose the function we're to investigate is a C(k+1) function within its open domain which includes a as an internal point. The purpose is to somehow find "the best" (polynomial) approximation to the function, at points like x, inside an open sphere around our reference point a, about which we know the value of the function there as well as all the (mixed) partial derivatives up to the kth degree. Now the first (and interstingly the most intuition-based) step to take is to somehow map our orthonormal bases to another orthonormal bases (using orthogonal transformations), such that one base will lie along the vector x-a. Now, with this change of coordinates, it's as though we're changing only one variable to go from a to x. Therefore, the problem turns into that of single variable Taylor expansion of the kth order. The other cool part of the attempt to "construct" the multivariable polynomial with only that for the single-variable at hand is that one will be able to and therefore should find all (mixed) partial derivatives of f at a, "this time in the direction of the vector x-a": (it becomes interesting when one thinks about how what they already have (here: all the (mixed) partials) enables them to obtain what they need to have (here: directional first, second,..., kth derivative of f at a with respect to x-a, if I'm rigth). It may also be illuminating to think about the relationship of the kth derivative in the direction of x-a and that of all the mixed partials at a.) Once the job is completed in thar coordinate system, one, quite clearly, needs to describe all they have obtained until there, in the language of the former coordinate system, thereby having constructed the multivariable polynomial approximation. 
A: You need to understand gradient vector to understand this article, although it is quite useful and can direct you to formulate the taylor series formula more concisely: http://www.math.ucdenver.edu/~esulliva/Calculus3/Taylor.pdf.
