Hello Abdeljalil!
Some comments, in line, below:



 abdeljalil:

I just got rid of torch.sigmoid() and used BCEWithLogitsLoss instead of BCELoss and it just worked fine as expected, there is no stuck loss problem anymore and the model converges within a couple of epochs.
…
But I’m not sure why this really solves my problem.


Mathematically, sigmoid() maps [-inf, inf] to [0.0, 1.0].
But sigmoid() has exp() in it, which “amplifies” numbers
dramatically, so numerically the range for which sigmoid()
gives useful results is rather small.
Specifically, sigmoid (x) saturates at 1.0 already for x = 17,
and saturates at 0.0 for x = -90 (for FloatTensors). These
are not extreme values and can easily be returned by a Linear
layer.
And once the sigmoid() saturates, its gradient becomes zero,
so the optimizer won’t know which direction to go to get you
“unstuck.”





 KFrank:

Note, when you convert to BCEWithLogitsLoss , when calculating
your accuracy, you should threshold the output of your model at
0.0 rather than 0.5 .


I really don’t know what you mean here exactly, but the model output now is just some integers, example output:
[ 7., -3., -5., -5., -5.,  4., -5., -5., -5.,  4., -5., 13., -5., -5., -5.,
-4.,  6.,  4.,  3.,  6., -3.,  3., -5.,  3., 11., -4., -4., 11., -5., -6.,
-5.,  4.,  6.,  7.,  5.,  7., -5.,  4.,  9., -5.,  8.,  4., 5.,  9., -4., 
1., 10.,  4., 14., -4., -4., -4.,  5., -5.,  6., 15., 3., -6., -5., -5.,
-5., -5., -4.,  3.]



These integer values for output surely occur only after you have
calculated output = torch.round (output).


So I just did the following:
output = torch.round(torch.sigmoid(output))



You can do this and it will give you what you want.  But, to me,
it’s a little more logically straightforward to compare (“threshold”)
your pre-sigmoid outputs with 0.
The rationale is:
An output of a linear layer (a raw-score logit) of 0.0 is mapped by
sigmoid() to a probability, P, of 0.5.  To convert a probability to
a discrete yes-no prediction,  one would typically say P > 0.5
means “yes” (and P < 0.5 means “no”).
So instead of using sigmoid() to convert to probabilities, you
can just threshold the logits directly: logit > 0.0 means “yes”
(and logit < 0.0 means “no”).


Yes running the same code in Keras always yields almost the same results, pretty robust.


Not that it really matters now that you have your model working,
but are you possibly using double-precision in your keras code?
Do you have an explicit sigmoid() in your keras code?  If so,
for what values does the keras sigmoid() saturate?
Best.
K. Frank