Avoiding subtraction for finite difference with log and exp I want to approximate the derivative of f(x)
Finite difference
$f'(x) \approx \frac{f(x+h)-f(x)}{h}$
I was taught that the error from the subtraction is blown up for small h. This I can verify with MATLAB.
Smartass method
So I though maybe the following would fix the problem:
$f'(x) \approx \frac{\log{\frac{\exp(f(x+h))}{\exp(f(x))}}}{h}$
However, I get the same relative errors. (Which also start increasing for h smaller than approx. $\epsilon_{mach}/2$
Why is this?
 A: Your solution doesn't help as it is essentially the same thing. $\log \frac{\exp(f(x+h))}{\exp(f(x))}$ is not going to be a more accurate representation of $f(x+h)-f(x)$ than evaluating it directly (on average, I would imagine it will be less accurate).
There isn't really much you can do when you are working that close to machine error afaik. Practically speaking, you don't need to use such a tiny $h$ with reasonable functions; the practical solution is that when you are dividing by a small number, don't let that number be close to $\epsilon_{mach}$.
A: First of all, you should always be careful to choose $h$ with a very slight dependence on $x$ such that $x \pm h$ is exactly machine representable. For instance, if $x = 1$, then you should at least round the $h$ you are tempted to use to the nearest multiple of $2^{-22}$ for 32-bit precision or $2^{-51}$ for 64-bit doubles.
Then, as Manazoni points out, the symmetric formula is superior; it gets a precision of about two orders of magnitude better working with 64-bit numbers. You should choose an $h$ somewhere near $x \epsilon^{1/3}$, where $\epsilon$ is the fractional machine precision
(thus in 64-bit doubles, $h$ of roughly $10^{-5}x$ is good.  The fractional error in the derivative should be about $\epsilon^{2/3}$.
Yet better, unless you suspect really weird behavior near $x$, is to use the next higher order symmetric formula
$$
f^{\prime}(x) = \frac{ f(x-2h) -4 f(x-h) + 4 f(x+h) - f(x+2h)}{6} + O\left(h^4 f{\prime\prime\prime\prime}(x) \right)
$$
using $h \approx x \epsilon^{1/5}$, which gives fractional accuracy of about $\epsilon^{4/5}$. 
Higher order methods begin to probe too sensitively for slightly irregular function behavior and begin to rely too heavily on accuracy of the function evaluations.  Unfortunately, unlike the case for integration where there are adaptive step size methods, there is no robust method of choosing order and step size that reliably gives you the best achievable 
estimate of the derivative. But there are methods that are yet better than the one I showed here, see chapter 5.7 of Numerical recipes.
