Estoy tratando de implementar el cálculo de gradiente para redes neuronales usando propagación hacia atrás. No puedo hacer que funcione con error de entropía cruzada y unidad lineal rectificada (ReLU) como activación.

Logré que mi implementación funcionara para el error al cuadrado con las funciones de activación sigmoid, tanh y ReLU. El error de entropía cruzada (CE) con gradiente de activación sigmoidea se calcula correctamente. Sin embargo, cuando cambio la activación a ReLU, falla. (Estoy omitiendo tanh para CE, ya que arroja valores en el rango (-1,1)).

¿Se debe al comportamiento de la función de registro en valores cercanos a 0 (que ReLU devuelve aproximadamente el 50% del tiempo para entradas normalizadas)? Traté de mitigar ese problema con:

log(max(y,eps))

pero solo ayudó a devolver el error y los gradientes a números reales; todavía son diferentes del gradiente numérico.

Verifico los resultados usando gradiente numérico:

num_grad = (f(W+epsilon) - f(W-epsilon)) / (2*epsilon)

El siguiente código de matlab presenta una implementación de retropropagación simplificada y condensada utilizada en mis experimentos:

function [f, df] = backprop(W, X, Y)
% W - weights
% X - input values
% Y - target values

act_type='relu';    % possible values: sigmoid / tanh / relu
error_type = 'CE';  % possible values: SE / CE

N=size(X,1); n_inp=size(X,2); n_hid=100; n_out=size(Y,2);
w1=reshape(W(1:n_hid*(n_inp+1)),n_hid,n_inp+1);
w2=reshape(W(n_hid*(n_inp+1)+1:end),n_out, n_hid+1);

% feedforward
X=[X ones(N,1)];
z2=X*w1'; a2=act(z2,act_type); a2=[a2 ones(N,1)];
z3=a2*w2'; y=act(z3,act_type);

if strcmp(error_type, 'CE')   % cross entropy error - logistic cost function
    f=-sum(sum( Y.*log(max(y,eps))+(1-Y).*log(max(1-y,eps)) ));
else % squared error
    f=0.5*sum(sum((y-Y).^2));
end

% backprop
if strcmp(error_type, 'CE')   % cross entropy error
    d3=y-Y;
else % squared error
    d3=(y-Y).*dact(z3,act_type);
end

df2=d3'*a2;
d2=d3*w2(:,1:end-1).*dact(z2,act_type);
df1=d2'*X;

df=[df1(:);df2(:)];

end

function f=act(z,type) % activation function
switch type
    case 'sigmoid'
        f=1./(1+exp(-z));
    case 'tanh'
        f=tanh(z);
    case 'relu'
        f=max(0,z);
end
end

function df=dact(z,type) % derivative of activation function
switch type
    case 'sigmoid'
        df=act(z,type).*(1-act(z,type));
    case 'tanh'
        df=1-act(z,type).^2;
    case 'relu'
        df=double(z>0);
end
end

Editar

Después de otra ronda de experimentos, descubrí que usar un softmax para la última capa:

y=bsxfun(@rdivide, exp(z3), sum(exp(z3),2));

y función de costo softmax:

f=-sum(sum(Y.*log(y)));

hacer que la implementación funcione para todas las funciones de activación, incluida ReLU.

Esto me lleva a la conclusión de que es la función de costo logístico (clasificador binario) que no funciona con ReLU:

f=-sum(sum( Y.*log(max(y,eps))+(1-Y).*log(max(1-y,eps)) ));

Sin embargo, todavía no puedo entender dónde radica el problema.