To numerically approximate the second partial derivatives, and get the Hessian matrix, one may use "2-point" central finite difference formulas:
$$H_{xy}(x,y) \approx \frac{f(x+h,y+h)-f(x+h,y-h)-f(x-h,y+h)+f(x-h,y-h)}{4h^2}\tag{1}$$
for the off-diagonal elements, and
$$H_{xx}(x,y) \approx \frac{f(x+h,y)-2f(x,y)+f(x-h,y)}{h^2}\tag{2}$$
formulas for the diagonal elements. (I am simplifying here, this generalizes for any number of dimensions, and uses the same displacement for all variables)
It may be sometimes desirable to further improve the accuracy of these central difference formulas. For univariate functions this is conveniently done via the 4-point formula (aka. 5-point stencil in one dimension):
$$f''(x) \approx \frac{-f(x+2 h)+16 f(x+h)-30 f(x) + 16 f(x-h) – f(x-2h)}{12 h^2}$$
However, I have been unable the find the "4-point" equivalent of $(1)$ and $(2)$ anywhere online.
Best Answer
For $H_{xx}$ it would be the same as the univariate $f''(x)$ you've posted, except for the extra argument for the other variable. $H_{yy}$ has the same form as $H_{xx}$ but with $x$ and $y$ swapped. The following assumes same step size for each dimension:
$$ \begin{aligned} H_{xx}(x,y) &\approx \frac{-f(x+2h,y)+16f(x+h,y)-30f(x,y)+16f(x-h,y)-f(x-2h,y)}{12h^2} \\ H_{yy}(x,y) &\approx \frac{-f(x,y+2h)+16f(x,y+h)-30f(x,y)+16f(x,y-h)-f(x,y-2h)}{12h^2} \\ \end{aligned} $$
For $H_{xy}$ it is
$$ H_{xy(x,y)} \approx \frac{ \begin{array}{c} -f(x+2h,y+2h)+16f(x+h,y+h)\\ +f(x+2h,y-2h)-16f(x+h,y-h)\\ +f(x-2h,y+2h)-16f(x-h,y+h)\\ -f(x-2h,y-2h)+16f(x-h,y-h)\\ \end{array} }{48h^2} $$
In case you're wondering where the formula comes from, you can derive it using Taylor series expansion.
$$ f(x+h)=f(x)+\frac{f'(x)}{1!}h+\frac{f''(x)}{2!}h^2 +\frac{f'''(x)}{3!}h^3+\cdots \tag{1}\label{eq1} $$
Substituting $h\to -h$ inverts the sign of each odd term in $\eqref{eq1}$
$$ f(x-h)=f(x)-\frac{f'(x)}{1!}h+\frac{f''(x)}{2!}h^2 -\frac{f'''(x)}{3!}h^3+\cdots \tag{2}\label{eq2} $$
Adding $\eqref{eq1}$ and $\eqref{eq2}$ together eliminates the odd terms
$$ f(x+h)+f(x-h)=2f(x)+2\frac{f''(x)}{2!}h^2 +2\frac{f^{(4)}(x)}{4!}h^4+\cdots \tag{3}\label{eq3} $$
Rearranging the terms, we arrive at
$$ \begin{aligned} f''(x)&=\frac{f(x+h)-2f(x)+f(x-h)}{h^2}-2\frac{f^{(4)}(x)}{4!}h^2+\cdots\\ &=\frac{f(x+h)-2f(x)+f(x-h)}{h^2}+O(h^2) \end{aligned} \tag{4}\label{eq4} $$
The terms from $f^{(4)}$ onwards are gathered into the truncation error term $O(h^2)$. This formulation is second-order accurate.
Substituting $h\to 2h$ in $\eqref{eq3}$ we get
$$ \begin{aligned} f(x+2h)+f(x-2h)&=2f(x)+2\frac{f''(x)}{2!}(2h)^2 +2\frac{f^{(4)}(x)}{4!}(2h)^4+\cdots\\ &=2f(x)+8\frac{f''(x)}{2!}h^2 +32\frac{f^{(4)}(x)}{4!}h^4+\cdots \end{aligned}\tag{5}\label{eq5} $$
Combining $\eqref{eq3}$ and $\eqref{eq5}$ with suitable coefficients we can eliminate $f^{(4)}$ term
$$ 16\eqref{eq3}-\eqref{eq5}=30f(x)+12h^2f''(x)-\frac{2}{15}h^6f^{(6)}(x)+\cdots \tag{6}\label{eq6} $$
which, after rearranging, we get a fourth-order accurate formulation for $f''$
$$ \begin{aligned} f''(x)&=\frac{-f(x+2h)+16f(x+h)-30f(x)+16f(x-h)-f(x-2h)}{12h^2}-\frac{2}{15}f^{(4)}(x)h^4+\cdots\\ &=\frac{-f(x+2h)+16f(x+h)-30f(x)+16f(x-h)-f(x-2h)}{12h^2}+O(h^4) \end{aligned}\tag{7}\label{eq7} $$
For bivariate function $f(x,y)$ the series expansion need to be applied separately for each variable. Otherwise, the line of thinking is the same: we expand the function on some points around $(x,y)$ and find a combination of coefficients for each point that will eliminate as much error terms as possible. The order of accuracy for the formula of $H_{xy}$ given in the question is $O(h^2)$, while the one in this answer is $O(h^4)$.
It is interesting to note that a different formulation exists for $H_{xy}$ which reuses $H_{xx}$ and $H_{yy}$
$$ H_{xy}=\frac{f(x+h,y+h)-2f(x,y)+f(x-h,y-h)}{2h^2}-\frac{1}{2}H_{xx}-\frac{1}{2}H_{yy}+O(h^2)\tag{8}\label{eq8} $$
Another one is $$ H_{xy}=\frac{-f(x+h,y-h)+2f(x,y)-f(x-h,y+h)}{2h^2}+\frac{1}{2}H_{xx}+\frac{1}{2}H_{yy}+O(h^2)\tag{9}\label{eq9} $$ while the 5-point stencil version is $$ H_{xy}=\frac{\begin{array}{l}-f(x+2h,y+2h)+16f(x+h,y+h)\\-30f(x,y)\\ +16f(x-h,y-h)-f(x-2h,y-2h)\end{array}}{24h^2}-\frac{1}{2}H_{xx}-\frac{1}{2}H_{yy}+O(h^4)\tag{10}\label{eq10} $$ and $$ H_{xy}=\frac{\begin{array}{l}f(x+2h,y+2h)-16f(x+h,y+h)\\+30f(x,y)\\ -16f(x-h,y-h)+f(x-2h,y-2h)\end{array}}{24h^2}+\frac{1}{2}H_{xx}+\frac{1}{2}H_{yy}+O(h^4)\tag{11}\label{eq11} $$ (in this case $H_{xx}$ and $H_{yy}$ need to be formulated using 5-point stencil as well) These formulations have the same order of accuracy, but the contents of the $O(h^2)$ and $O(h^4)$ are different, and may have slightly larger error. The benefit of using this alternative formulation is the reduced number of points required to calculate the Hessian (9 points reduced to 7 for second order, 17 points reduced to 13 points for fourth order accuracy).