Choosing cost function for a neural network continuous forecast I´m using Octave to build a Neural Network regression (3 layers nn, using tanh as activation function, fmincg as optimizer and continuous output). The purpose of the model is to forecast demand for some products, based on multiple variables (past demand of the products, past stock levels, pendent orders, etc.). I'm normalizing all inputs and output.
The quality of the forecast is finally evalued with SFE (statistical forecast error), which I defined in Octave as:
    $$SFE=\frac{\sum_{i=1}^{m}abs(Y_m-Forecast_m)}{\sum_{i=1}^{m}Forecast_m}$$.
Despite this, I'm training the model as a "traditional" regression, using a sum of square errors as cost function (J):
    $$J=\frac{\sum_{i=1}^{m}(Y_m-Forecast_m)^2}{2m}$$.
I think this is not correct, as learning should be made towards the real objective (SFE) and J is over-estimating the impact of bigger deviations. But as abs() has no derivative, I don't realize how SFE function could work under a back-prop training for the nn.
How should the cost function (J) be defined so as to train the model towards SFE? 
I'd like to use SFE itself, but I don't find it possible as I have to compute the gradient in order for the optimizer to work (fmincg for Octave) and abs() is not differentiable at 0. The answer should cover this aspect for the optimizer to work.
 A: It would probably be best to use a loss function based on the absolute, rather than the squared errors.  This will give a model that computes the conditional median of the target distribution, rather than the conditional mean.  For details, see
M Saerens, "Building cost functions minimizing to some summary statistics", IEEE Transactions on Neural Networks, 11 (6), 1263-1271, 2000.
which is a very nice paper and ought to be much better known that it is.  It is probably not too doficult to use the SFE itself as the loss function, although it will involve more programming.
A: The $L^2$ (or quadratic) error is defined as
$$J_2 = \frac{1}{m} \sum_{i=1}^m {(Y_i - F_i)^2}$$
whereas the $L^1$ (or absolute) error is defined as
$$J_1 = \frac{1}{m} \sum_{i=1}^m{ |Y_i - F_i|}$$
$J_1$ should have similar properties to $SFE$, though the gradients w.r. to $F_i$ will be easier to compute as there is no forecast dependent denominator.
The gradient of the $L^1$ norm $|x|$ w.r. to $x$ can be taken to be 1 for positive values and -1 for negative values. You will find that the special case $x=0$ does really not matter numerically, as it is unlikely you hit 0, so you can use -1, 1 or 0 for this value.
This suggests:
$$\partial_i J_1 = \frac{1}{m} (2 \cdot 1_{\{Y_i>F_i\}} - 1)$$
