Why cant I overfit this dataset with my neural network? I have read that given a model is complex enough and I train for enough epochs, my model should at some point overfit the dataset. However I implemented a simple neural network in keras and my validation set loss seems to never go up:
import tensorflow as tf
from tensorflow import keras
import numpy as np
import random
from sklearn import datasets, preprocessing
import matplotlib.pyplot as plt


# import and scale
dataset = datasets.load_boston()
X = dataset.data
y = dataset.target
X = preprocessing.scale(X)
y = y.reshape((y.shape[0], 1))

# shuffle
shuffle_indices = list(range(X.shape[0]))
random.shuffle(shuffle_indices)
X = X[shuffle_indices]
y = y[shuffle_indices]

# tain-validation split
X_train, y_train  = X[:int(X.shape[0] * 0.7)], y[:int(X.shape[0] * 0.7)]
X_val, y_val = X[int(X.shape[0] * 0.7):], y[int(X.shape[0] * 0.7):]

# define and fit model
model = keras.Sequential([keras.layers.Dense(X.shape[1], use_bias=True, activation=&quot;sigmoid&quot;),
                              keras.layers.Dense(128, use_bias=True, activation=&quot;sigmoid&quot;),
                              keras.layers.Dense(128, use_bias=True, activation=&quot;sigmoid&quot;),
                              keras.layers.Dense(128, use_bias=True, activation=&quot;sigmoid&quot;),
                              keras.layers.Dense(128, use_bias=True, activation=&quot;sigmoid&quot;),
                              keras.layers.Dense(128, use_bias=True, activation=&quot;sigmoid&quot;),
                              keras.layers.Dense(128, use_bias=True, activation=&quot;sigmoid&quot;),
                              keras.layers.Dense(y.shape[1])
                         ])
model.compile(optimizer=tf.keras.optimizers.SGD(
    learning_rate=0.0001
), loss='MeanSquaredError')

model.fit(X_train, y_train, validation_data=(X_val, y_val), epochs=500, verbose=1)

# plot data
plt.plot(range(1, len(model.history.history['loss']) + 1), model.history.history['loss'], label='Train Set Cost')
plt.plot(range(1, len(model.history.history['val_loss']) + 1), model.history.history['val_loss'], label='Validation Set Cost')
plt.xlabel(&quot;epoch&quot;)
plt.ylabel(&quot;loss&quot;)
plt.legend()
plt.show()

The model is a simple dense neural network with Mean Squared Error as its loss function and gradient descent as it's optimizer. I tried to make the network deeper, but the validation loss only keeps decreasing until it stops at one point.
It is usually good to start with a small model because you can then evaluate the contribution of adding layers, etc.
Also, Boston dataset is a popular dataset so there are several tutorials showing good neural network architectures, like this one.
Concerning your model, here are some notes.

The use of sigmoid activation is likely to worsen results, since sigmoid function compresses values between 0 and 1, while you are trying to predict outputs between 5 and 50.
Instead of sigmoid, you can use ReLU activation, which has better convergence properties for inner layers
You can standardize your target data to reduce the variance of your data and control the mean, this usually improve a lot regression models. In this case, a sigmoid activation on your final layer would be a good choice