function subgrad = constraint_subgradient(Ac, idmalec, idfemalec, nummalec, numfemalec, w)

    % compute mean_male = mean( sigmoid ( A_male * w ) )
    sigmoid_mean_male = mean( 1 ./ ( 1 + exp( - Ac(idmalec,:)*w ) ) );
    % compute mean_female = mean( sigmoid ( A_female * w ) )
    sigmoid_mean_female = mean( 1 ./ ( 1 + exp( - Ac(idfemalec,:)*w ) ) );
    % compute the first part, sgn value = sgn ( abs( mean_male - mean_female ) ), in the subgradient G'(w)
    sgn = sign( sigmoid_mean_male - sigmoid_mean_female );

    % compute the second part in the subgradient G'(w)
    % compute derivative_male = Ac_male^T * ( sigmoid ./ (1 - sigmoid ) ) / nc_M
    sigmoid_male = 1 ./ ( 1 + exp( - Ac(idmalec,:)*w ) );
    subgrad_male = Ac(idmalec,:)' * ( sigmoid_male .* (1 - sigmoid_male) ) / nummalec;
    % compute derivative_female = Ac_female^T * ( sigmoid ./ (1 - sigmoid ) ) ) / nc_F
    sigmoid_female = 1 ./ ( 1 + exp( - Ac(idfemalec,:)*w ) );
    subgrad_female = Ac(idfemalec,:)' * ( sigmoid_female .* (1 - sigmoid_female) ) / numfemalec;

    % get the value of subgradient G'(w)
    subgrad = sgn * ( subgrad_male - subgrad_female );
end