Pytorch: NN function approximator, 2 in 1 out
[Please be aware of the Edit History below, as the major problem statement has changed.]

We are trying to implement a neural network in pytorch, that approximates a function f(x,y)=z. So there are two real numbers as input and one as ouput, we therefore want 2 nodes in the input layer and one in the output layer. We constructed a test set of 5050 samples and had pretty good results for that task in Keras with Tensorflow backend, with 3 hidden layers with a configuration of the nodes like: 2(in) - 4 - 16 - 4 - 1(out); and ReLU activation functions on all hidden layers, linear on in- and output.

Now in Pytorch we tried to implement a similar network but our loss function still literally explodes: It changes in the first few steps and converges then to some value around 10^7. In Keras we had an error around 10 percent. We already tried different network configurations without any improvement. Maybe someone could have a look on our code and suggest any change?

To explain: tr_data is a list, containing 5050 2*1 numpy arrays which are the inputs for the network. tr_labels is a list, containing 5050 numbers which are the outputs we want to learn. loadData() just load those two lists.

import torch.nn as nn
import torch.nn.functional as F

BATCH_SIZE     = 5050
DIM_IN         = 2
DIM_HIDDEN_1   = 4
DIM_HIDDEN_2   = 16
DIM_HIDDEN_3   = 4
DIM_OUT        = 1
LEARN_RATE     = 1e-4
EPOCH_NUM      = 500

class Net(nn.Module):
    def __init__(self):
        #super(Net, self).__init__()
        super().__init__()
        self.hidden1 = nn.Linear(DIM_IN, DIM_HIDDEN_1)
        self.hidden2 = nn.Linear(DIM_HIDDEN_1, DIM_HIDDEN_2)
        self.hidden3 = nn.Linear(DIM_HIDDEN_2, DIM_HIDDEN_3)
        self.out     = nn.Linear(DIM_HIDDEN_3, DIM_OUT)

    def forward(self, x):
        x = F.relu(self.hidden1(x))
        x = F.tanh(self.hidden2(x))
        x = F.tanh(self.hidden3(x))
        x = self.out(x)
        return x

model = Net()

loss_fn = nn.MSELoss(size_average=False)

optimizer = torch.optim.Adam(model.parameters(), lr=LEARN_RATE)

tr_data,tr_labels = loadData()

tr_data_torch    = torch.zeros(BATCH_SIZE, DIM_IN)
tr_labels_torch  = torch.zeros(BATCH_SIZE, DIM_OUT)

for i in range(BATCH_SIZE):
    tr_data_torch[i]   = torch.from_numpy(tr_data[i])
    tr_labels_torch[i] = tr_labels[i]

for t in range(EPOCH_NUM):

    labels_pred = model(tr_data_torch)
    loss = loss_fn(labels_pred, tr_labels_torch)
    #print(t, loss.item())

    optimizer.zero_grad()
    loss.backward()
    optimizer.step()


I have to say, those are our first steps in Pytorch, so please forgive me if there are some obvious, dumb mistakes. I appreciate any help or hint,
Thank you! 

EDIT 1 ------------------------------------------------------------------

Following the comments and answers, we improved our code. The Loss function has now for the first time reasonable values, around 250. Our new class definition looks like:

class Net(nn.Module):
    def __init__(self):
        super(Net, self).__init__()
        #super().__init__()
        self.hidden1 = nn.Sequential(nn.Linear(DIM_IN, DIM_HIDDEN_1), nn.ReLU())
        self.hidden2 = nn.Sequential(nn.Linear(DIM_HIDDEN_1, DIM_HIDDEN_2), nn.ReLU())
        self.hidden3 = nn.Sequential(nn.Linear(DIM_HIDDEN_2, DIM_HIDDEN_3), nn.ReLU())
        self.out = nn.Linear(DIM_HIDDEN_3, DIM_OUT)

    def forward(self, x):
        x = self.hidden1(x)
        x = self.hidden2(x)
        x = self.hidden3(x)
        x = self.out(x)
        return x


and the loss function:

loss_fn = nn.MSELoss(size_average=True, reduce=True)


As we stated before, we already had far more satisfying results in keras with tensorflow backend. The loss function was around 30, with a similar network configuration. I share the essential parts(!) of our keras code here:

model = Sequential()
model.add(Dense(4, activation="linear", input_shape=(2,)))
model.add(Dense(16, activation="relu"))
model.add(Dense(4, activation="relu"))
model.add(Dense(1, activation="linear" ))
model.summary()

model.compile ( loss="mean_squared_error", optimizer="adam", metrics=["mse"] )

history=model.fit ( np.array(tr_data), np.array(tr_labels), \
                    validation_data = ( np.array(val_data), np.array(val_labels) ), 
                    batch_size=50, epochs=200, callbacks = [ cbk ] )


Thank your already for all the help! If anybody still has suggestions to improve the network, we would be happy about it. As somebody already asked for the data, we want to share a pickle file here: 

https://mega.nz/#!RDYxSYLY!P4a9mEDtZ7A5Bl7ZRjRk8EzLXQt2gyURa3wN3NCWFPA

together with the code to access it:

import pickle
f=open("data.pcl","rb")
tr_data=pickle.load ( f )
tr_labels=pickle.load ( f )
val_data=pickle.load ( f )
val_labels=pickle.load ( f )
f.close()


It should be interesting for you to point out the differences between torch.nn and torch.nn.functional (see here). Essentially, it might be that your backpropagation graph might be executed not 100% correct due to a different specification.

As pointed out by previous commenters, I would suggest to define your layers including the activations. My personal favorite way is to use nn.Sequential(), which allows you to specify multiple opeations chained together, like so:

self.hidden1 = nn.Sequential(nn.Linear(DIM_IN, DIM_HIDDEN1), nn.ReLU())


and then simply calling self.hidden1 later (without wrapping it in F.relu()). 

May I also ask why you do not call the commented super(Net, self).__init__() (which is the generally recommended way)?
Additionally, if that should not fix the problem, can you maybe just share the code for Keras in comparison?
