As a solution, you can train single layer linear perceptron to find weights of each feature. After that iteratively drop the most useless feature (with min weight) and test combinations of remained features (with the highest accuracy) on your NN. Code below
def features_looker(X_tr, y_tr, X_ts, y_ts, remove_features=0):
model = tf.keras.Sequential()
model.add(tf.keras.layers.Dense(units=1, input_dim=X_tr.shape[1], activation='sigmoid'))
model.compile(loss='mean_squared_error', optimizer='adam', metrics=['accuracy'])
info = model.fit(X_tr, y_tr, epochs=30, batch_size=100, verbose=0, validation_data=(X_ts, y_ts))
accuracy = model.evaluate(X_tr, y_tr, verbose=0)[1]
val_accuracy = model.evaluate(X_ts, y_ts, verbose=0)[1]
col = X_tr.columns
weights = np.abs(model.layers[0].get_weights()[0].flatten())
weights_df = pd.DataFrame({'feature': col, 'weight': weights}).sort_values(by='weight', ascending=False)
worst_features = weights_df.nsmallest(remove_features, 'weight')['feature'].values
print('Features: ', len(col))
print('Train Accur: ', accuracy)
print('Test Accur: ', val_accuracy)
return info, worst_features, weights_df