Taylor Expansion for a Multivariable Function \begin{align}
f(x_1,\dots,x_d) &= \sum_{n_1=0}^\infty \sum_{n_2=0}^\infty \cdots \sum_{n_d = 0}^\infty 
\frac{(x_1-a_1)^{n_1}\cdots (x_d-a_d)^{n_d}}{n_1!\cdots n_d!}\,\left(\frac{\partial^{n_1 + \cdots + n_d}f}{\partial x_1^{n_1}\cdots \partial x_d^{n_d}}\right)(a_1,\dots,a_d) \\
&= f(a_1, \dots,a_d) \\
&\quad + \sum_{j=1}^d \frac{\partial f(a_1, \dots,a_d)}{\partial x_j} (x_j - a_j) \\ 
&\quad + \sum_{j=1}^d \sum_{k=1}^d \frac{1}{2!} \frac{\partial^2 f(a_1, \dots,a_d)}{\partial x_j \partial x_k} (x_j - a_j)(x_k - a_k) \\ 
&\quad + \sum_{j=1}^d\sum_{k=1}^d\sum_{l=1}^d \frac{1}{3!} \frac{\partial^3 f(a_1, \dots,a_d)}{\partial x_j \partial x_k \partial x_l} (x_j - a_j)(x_k - a_k)(x_l - a_l) \\ 
&\quad + \dots
\end{align}
I have read through wikipedia, and when I saw the above formula, I didn't know how the second equality is justified. Anybody can help me please? (The first equality is assumed to be true by myself thus doesn't need to be proved)
 A: The second RHS is an enumeration of the first RHS according to the value of $m=n_1+\cdots+n_d$. For $m=0$, one gets one term, which is $f(a_1, \dots,a_d)$. For $m=1$, one gets $d$ terms, which are the products $\frac{\partial f(a_1, \dots,a_d)}{\partial x_j}\cdot(x_j - a_j)$ for each $1\leqslant j\leqslant d$.
More generally, for each $m\geqslant0$, one gets $d^m$ terms, hence the multiple sums from $1$ to $d$ with $m$ sums.
To "sum" the above, one uses the identity
$$
\sum_{n_1=0}^\infty \sum_{n_2=0}^\infty \cdots \sum_{n_d = 0}^\infty A(n_1,\cdots,n_d) =\sum_{m=0}^\infty\sum_{\begin{array}{c}(n_1,\cdots,n_d)\\ n_1+\cdots+n_d=m\end{array}} A(n_1,\cdots,n_d),
$$
with
$$
A(n_1,\cdots,n_d)=\frac{\partial^m f(a_1, \dots,a_d)}{\partial^{n_1} x_1\cdots \partial^{n_d} x_d}\cdot\prod_{j=1}^m (x_j - a_j)^{n_j}.
$$
A: Let me just take an analogy. When you make a first order Taylor expansion of $f(x)$, you basically write the equation of a straight line saying that the model is linear with respect to $x$. If you do the same with $g(x,y)$ and you want the model to be linear with respect to both $x$ and $y$, you need to write that $$g(x,y) \simeq a +b (x-x_0)+c(y-y_0)+d(x-x_0)(y-y_0)$$  
If you prefer, say that for a given value of $y$,$\text{  } g(x,y)$ is linear with respect to $x$; this write $$g(x,y)=a(y)+b(y) \times (x-x_0)$$ and now consider that $a(y)$ and $b(y)$ are expanded as Taylor series around $y_0$. So $$a(y)=\alpha_0+\alpha_1 (y-y_0)$$ $$b(y)=\beta_0+\beta_1 (y-y_0)$$ and replace in the previous expansion for $g(x,y)$.
You can generalize this for as many variables as you wish. I let you finding the analogy between the coefficients and the derivatives.
A: For simplicity, I demonstrate the two-dimensional case, up to the second order. The multi-dimensional case is similar.
Suppose that the function $f=f(x,y)$ is defined on a convex open
subset of $\mathbb{R}^{2}$, and suppose $f$ is sufficiently smooth. Let $(x_{0},y_{0})$
be a point in the domain of $f$. Let $h,k\in\mathbb{R}$ such that $(x_{0}+h,y_{0}+k)$
is also point in the domain of $f$. Define a function $g(t)=f(x_{0}+th,y_{0}+tk)$.
Note that $g$ is well-defined on an open interval containing $[0,1]$
and $g$ is sufficiently differentiable. Consider the Taylor expansion
of the one-variable function $g$: 
$$
f(x_{0}+h,y_{0}+k)-f(x_{0},y_{0})=g(1)-g(0)=g'(0)(1-0)+\frac{1}{2!}g''(0)(1-0)^{2}+\ldots
$$
By chain rule, $g'(t)=f_{x}h+f_{y}k,$ $g''(t)=(f_{xx}h+f_{xy}k)h+(f_{yx}h+f_{yy}k)k=f_{xx}h^{2}+2f_{xy}hk+f_{yy}k^{2}$,
where the derivatives $f_{x},f_{y}, f_{xx}$ etc are evaluated at the point
$(x_{0}+th,y_{0}+tk)$. Plug-in $t=0$, we have $g'(0)=f_{x}(x_{0},y_{0})h+f_{y}(x_{0},y_{0})k$
and $g''(0)=\ldots$ (I don't want to type). The result follows. 
