Second derivative of a composition? Let $g:\mathbb{R}^n\to\mathbb{R}^p$, $f:\mathbb{R}^p\to\mathbb{R}$ and define the composition $h(x) = f(g(x))$. The gradient of $h$ with respect to $x$, $\nabla h\in\mathbb{R}^{1\times n}$, is given by
\begin{align}
\nabla h(x) = \underbrace{\nabla f(g(x))}_{1\times p} \underbrace{\nabla g(x)}_{p\times n}
\end{align}
How do I compute $\nabla^2 h(x)$?
 A: Clearly the first derivative is $D(f\circ g)_x;h\in\mathbb{R}^n\rightarrow Df_{g(x)}(Dg_x(h))\in\mathbb{R}$. Then the second derivative is
$D^2(f\circ g)_x:(h,k)\in(\mathbb{R}^n)^2\rightarrow D^2f_{g(x)}(Dg_x(h),Dg_x(k))+Df_{g(x)}(D^2g_x(h,k))\in\mathbb{R}$.
Now, if we want to use matricial notations: $[D^2(f\circ g)_x]$ is a $(n\times n)$ symmetric matrix, $[D^2f_{g(x)}]$ is a symmetric $(p\times p)$ matrix, $[Dg_x(h)],[Dg_x(k)]\in M_{p1}$, $[D^2g_x(h,k)]\in M_{p1}$ and $[Df_{g(x)}]=\nabla f(g(x))\in M_{1,p}$. Moreover $[Dg_x]\in M_{pn}$.
We obtain $h^T[D^2(f\circ g)_x]k=h^T[Dg_x]^T[D^2f_{g(x)}][Dg_x]k+\nabla f(g(x))[D^2g_x(h,k)]$.
EDIT 1. Of course, $[D^2g_x]$ is more complicated: it is a stack of $p$ $(n\times n)$ symmetric matrices: $[D^2g_x]=[D^2g_{1x},\cdots,D^2g_{px}]^T$ and $[D^2g_x(h,k)]=[h^TD^2g_{1x}k,\cdots,h^TD^2g_{px}k]^T$. Thus $\nabla f(g(x))[D^2g_x(h,k)]=\sum_{i=1}^p \dfrac{\partial f}{\partial y_i}(g(x))h^TD^2g_{ix}k=h^T(\sum_{i=1}^p \dfrac{\partial f}{\partial y_i}(g(x))D^2g_{ix})k=h^T(\nabla f(g(x))[D^2g_x])k$.
Conclusion: $[D^2(f\circ g)_x]=[Dg_x]^T[D^2f_{g(x)}][Dg_x]+\nabla f(g(x))[D^2g_x]$. Note that the above matrix is symmetric.
EDIT 2. Answer to @ rogerG . 1. The first derivative is a linear application: $Dg_x\in L(\mathbb{R}^n,\mathbb{R}^p),Df_{g(x)}\in L(\mathbb{R}^p,\mathbb{R}),D(f\circ g)_x=Df_{g(x)}\circ Dg_x\in L(\mathbb{R}^n,\mathbb{R})$.


*The second derivative of a function $\phi:\mathbb{R}^n\rightarrow \mathbb{R}$ is a symmetric bilinear application: $D^2\phi_x\in BL(\mathbb{R}^n\times\mathbb{R}^n\rightarrow \mathbb{R})$. Its associated matrix is s.t. $[D^2\phi_x]_{ij}=\dfrac{\partial^2\phi}{\partial x_ix_j}(x)$.

*Your required formula is a generalization of the formula when $n=p=1$: $(f\circ g)'=(f'\circ g)g'$ and $(f\circ g)"=(f"\circ g)g'^2+(f'\circ g)g"$.
