Which is the correct Hessian matrix (the standard matrix of a bilinear form)? Please note the typo in the first entry: $\frac {\partial^2f} {\partial x_1\partial x_2}$ should instead be $\frac{\partial^2f} {\partial x_1\partial x_1}$.
Also, this Hessian matrix need not be symmetric as the partials need not all be continuous.
The convention in use is: $\frac {\partial^2f} {\partial x_1\partial x_2}$ means $\partial_{x_1} (\partial_{x_2}f)$.

Question: Shouldn't the correct matrix be the transpose of what is
given?
My reasoning is this:
Let $\{e_1, ..., e_n\}$ be the standard basis for $\mathbb R^n$.
Because the matrix is for the bilinear form
$D^2f(x)$, by definition, the $ij$-entry of the matrix is given by
$$[D^2f(x)](e_i,e_j)=[[D(Df)(x)] (e_i)](e_j) =\frac {\partial^2f} {\partial x_j\partial x_i} \neq \frac {\partial^2f} {\partial x_i\partial x_j}.$$ In other words, in the $ij$-entry, the $i$th variable should be partialled before the $j$th variable, not the other way round.

Referring to the proof is not helpful as it is cryptic and makes no sense:


Information on how $D^2f(x)$ is defined:
$\newcommand{\hom}{\operatorname{Hom}}$Let us write $A \subseteq U=V=\mathbb{R}^n$. We have
\begin{align*}
f:&amp; A \to \mathbb{R}\\
Df:&amp; A \to \hom(V,\mathbb{R})\\
D^2f:=D(Df):&amp; A \to \hom(U, \hom(V,\mathbb{R}))
\end{align*}
Note that the codomain of $Df$ is $\hom(V,\mathbb{R})$ (the set of all linear maps from $V$ to $\mathbb{R}$), which is not a Euclidean space. So, what do we mean by $D(Df)$? Since $\hom(V,\mathbb{R})$ is isomorphic to $\mathbb{R}^n$, the answer is that we first identify each linear operator $L_i: e_j\mapsto\delta_{ij}$ with the vector $e_i$. Hence the linear map $Df(x)$ is identified with the gradient vector $\nabla f(x)=\left(\frac{\partial f}{\partial x_1},\,\frac{\partial f}{\partial x_2},\,\ldots,\,\frac{\partial f}{\partial x_n}\right)^\top$ before taking a second derivative.

Having this identification, the matrix of $D^2f$, as a linear operator, is the Jacobian matrix of $\nabla f$, i.e. $J(\nabla f)(x)=\left(\frac{\partial^2 f}{\partial x_j \partial x_i}\right)$. Therefore, $D^2f(x)(u)$ is represented by the vector $\left(\frac{\partial^2 f}{\partial x_j \partial x_i}\right)u$ and
$$
\left(D^2f(x)(u)\right)(v) = \left[\left(\frac{\partial^2 f}{\partial x_j \partial x_i}\right)u\right]^\top v = u^\top \left(\frac{\partial^2 f}{\partial x_i \partial x_j}\right)v.\tag{1}
$$
Technically, $D^2f(x)$ is not a bilinear form, but a linear operator that maps vectors to linear operators. Yet, the mapping $B_x:U\times V\to\mathbb{R}$ by $B_x(u,v)=\left(D^2f(x)(u)\right)(v)$ is bilinear. Therefore we can identify $D^2f(x)$ with the bilinear form $B_x$. And when we speak of the matrix of $D^2f(x):U\times V\to\mathbb{R}$, we actually mean the matrix of $B_x$ (and $D^2f(x)$ is not really a function from $U\times V$ to $\mathbb{R}$ in the first place).

By convention, if $\mathbf{u},\mathbf{v}$ and $\mathbf{B}$ represent respectively two vectors $u,v$ and a bilinear form $b(u,v)$ with respect to some basis, then $b(u,v)=\mathbf{u}^\top \mathbf{B}\mathbf{v}$. Therefore, from $(1)$, we see that the matrix of $B_x$ with respect to the standard basis is $\left(\frac{\partial^2 f}{\partial x_i \partial x_j}\right)$.
If $f$ is twice differentiable (in sense of totally / Frechet differentiable), then the Hessian is symmetric.

Notice that there is little merit to using the Hessian if $f$ is not twice differentiable but only twice partially twice differentiable.

Proof of Schwarz's theorem (for twice differentiable function)

Without loss of generality assume the domain $A = [-1, 1]^2$ is a square, $f$ differentiable on $A$, $\nabla f$ is differentiable at $0$, and $f(0) = 0$ and $\nabla f(0) = 0$. We show that $\partial_y \partial_x f(0) = \partial_x \partial_y f(0)$.
The proof is little bit ugly. The key idea is that you can expand $f(x,y)$ in two way:
\begin{align} 
f(x,y) &amp;= \frac12 \partial_x \partial_x f(0) x^2 + \frac12 \partial_y \partial_y f(0) y^2 +  \partial_y \partial_x f(0) x y + \mathcal o(\|(x,y)\|^2), \\
f(x,y) &amp;= \frac12 \partial_x \partial_x f(0) x^2 + \frac12 \partial_y \partial_y f(0) y^2 +  \partial_x \partial_y f(0) x y + \mathcal o(\|(x,y)\|^2).
\end{align}
Thus, it follows
$$ \partial_y \partial_x f(0) - \partial_x \partial_y f(0) = \mathcal o(1). $$
That is, $\partial_y \partial_x f(0) = \partial_x \partial_y f(0)$.

Let us see, how to get the first expansion. 
Since $\nabla f$ is continuous at $0$, we may assume $\nabla f$ is bounded on $A$ without loss of generality.
Thus, by the Taylor theorem, we have
\begin{align}
f(x,y) &amp;= f(0,y) + \int_0^1 \partial_x f(t x, y) \;\mathrm dt \cdot x \\
&amp;= \frac12 \partial_y\partial_y f(0) y^2 + \mathcal o(y^2) + \int_0^1 \partial_x f(t x, y) \;\mathrm dt \cdot x \\
&amp;= \frac12 \partial_y\partial_y f(0) y^2 + \mathcal o(y^2) + \int_0^1 \partial_x \partial_x f(0)tx + \partial_y \partial_x f(0)y \;\mathrm dt \cdot x + \mathcal o(x^2 + y^2) \\
&amp;= \frac12 \partial_y\partial_y f(0) y^2 + \frac12\partial_x \partial_x f(0)x^2 + \partial_y \partial_x f(0)xy + \mathcal o(x^2 + y^2).
\end{align}
Analogously, we have
\begin{align}
f(x,y) &amp;= f(x,0) + \int_0^1 \partial_y f(x, ty) \;\mathrm dt \cdot y \\
&amp;= \frac12 \partial_x\partial_x f(0) x^2 + \frac12\partial_y \partial_y f(0)y^2 + \partial_x \partial_y f(0)xy + \mathcal o(x^2 + y^2).
\end{align}
Which had to be proven.