Why is central difference preferred over backward and forward difference in convolution? It is mentioned in some literature that we should always use central difference when computing the derivatives of an image instead of forward or backward difference. Does anyone knows why is that? 
Central difference = $\frac{df(x)}{dx} = \frac{f(x+h) - f(x-h)}{2h}$ 
Forward difference = $\frac{df(x)}{dx} = \frac{f(x+h) - f(x)}{h}$ 
Backward difference = $\frac{df(x)}{dx} = \frac{f(x) - f(x-h)}{h}$ 
 A: A few things come to mind:


*

*For smooth $f$, the central difference scheme is second order in $h$, whereas the other two you mentioned are first order in $h$. In other words, if $f$ is smooth, the (real space) error for the centered difference scheme is $O(h^2)$ whereas for the forward/backward schemes it is $O(h)$. Thus at fixed sufficiently small $h$, the central difference will have a better (real space) error than the other two. (Specifically, if $f \leq c_1 h$ and $g \leq c_2 h^2$ then $g<f$ once $h<\frac{c_2}{c_1}$. In this context, $c_1$ is proportional to a bound on $f''$ while $c_2$ is proportional to a bound on $f'''$.)

*Central differencing uses the same number of points as the other two you mentioned, so there is no loss in efficiency compared to those.

*There are higher order methods even than central differencing, and actually some of these may be even better. But very high order methods are not practical in floating point arithmetic, because they get their high order from subtraction of large nearly equal numbers, which causes loss of precision. "Middle order" methods (of order 2 through 6) are generally preferred when using double precision, because they do a good job of balancing accuracy in exact arithmetic with roundoff error when the calculations are done in double precision. 


There are some other aspects of this that are more specific to signal processing. If this doesn't answer your question by itself then I can mention a little bit about that.
Upon request, here is a figure comparing the (signed) error in the central, forward, and backward derivative estimates for $e^x$ at $x=1$. The centered difference is in black, the forward difference is in blue, and the backward difference is in red.

A: You already have got a couple of good relevant points, so I'm just gonna add one I haven't seen so far among the answers.


*

*The result of an operator with a well defined center pixel is on the same grid where you could argue that forward or backward difference are off by a fraction of 1/2 samples in either dimension (compared to the in-grid), this  could be impractical for many reasons and in many circumstances.


Here we can see how the grid moves from before and to after forward x-difference and forward y-difference respectively: None of the 3 grids are the same!

And for the central difference scheme all of them align nicely on top of each other: 
And as many things in engineering and science estimating a differential is one of many operations which need to be compatible or align in some sense for them to be useful.

EDIT
As proposed by @Ian in the comment to avoid the grids not being aligned, one could use (expressed in the functions Z-transform):
$$\mathcal{Z}\{d_k(x_1,x_2)\} = 1-{z_k}^2$$
But we see that it will simply be a lazy filtering by one grid step of the central difference:
$$\mathcal{Z}\{d_k(x_1,x_2)\} = \underset{\text{lazy}}{\underbrace{z_k}}\underset{\text{central diff.}}{\underbrace{({z_k}^{-1}-{z_k}^1)}}$$
So in practice we would be calculating the exact same thing, just moving the result one step to the left (or down).
A: One of the reasons, not unique to image processing, is that the central difference produces a 2nd order error, while the others - 1st order, so usually the central difference schemes have better convergence and/or stability properties.
