I don't understand what the fuzz is about all this forward-mode AD stuff. This is the easy part. Hell, you can basically use complex numbers to do all this for you for free (use Cauchy-rule + FD) without changing the rest of your code.
If that is what they meant, that's not a good idea.
Choice of h dictates precision (less bias than doing it purely in the real domain, but still .. a nontrivial choice). And also requires your functions to be analytical - which many aren't.
I know the contribution of the higher order imaginary terms are negligible in the result, but dual numbers are the more appropriate choice in this instance correct? since the higher order terms will vanish.