Why doesnt my lstm model for time series prediction improve after certain level of performance?

Question

I created an lstm model which predicts multioutput sequeances. It takes variable length sequences as input. These sequences are padded with zero to obtain equal length. Note that the time series are not equally spaced but time stamp is added as predictor. The descrete time series predictors are normalized with mean and standard deviation and run through PCA, the categorical feature is one hot encoded and the ordinal features are integer encoded. So the feature set is a combination of dynamic and static features. The targets are also scaled between [-1,1]. The input layer goes through a masked layer and then to 2 stacked lstm layers and dense layer to predict the targets (see code below).

Training actually starts well but then the performance starts to saturate. It seems the network also focus more on the 3rd output rather than the first two. This is seen in the validation curve of the third output that follows the training curve perfectly. For the first 2 outputs, the network has a hard time predicting some peak values. I have been tuning the hyper parameters but validation error does not go below a certain value. The longer I train the more the validation curve and training curve separate from each other and overfitting occurs on the first 2 outputs. I tried all the standard initializations and he_initialization seems to work the best. When more data is added there is a slight improvement in validation error but not significant. When adding dropout the validation error is lower than the training error due to noise introduced by dropout in feed forward but there is no significant improvement. Since neural networks tend to converge close to where they are initialized, I was thinking my initialization is not good.

I was wandering if anyone had any suggestions on how to improve the error of this model. I think I will be happy if I can get the validation error somewhere around 0.01.

def masked_mse(y_true, y_pred):
    mask = keras.backend.all(keras.backend.not_equal(y_true, 0.), axis=-1, keepdims=True)

    y_true_ = tf.boolean_mask(y_true, mask)
    y_pred_ = tf.boolean_mask(y_pred, mask)

    return keras.backend.mean(keras.backend.square(y_pred_ - y_true_))

def rmse(y_true, y_pred):
    # find timesteps where mask values is not 0.0
    mask = keras.backend.all(keras.backend.not_equal(y_true, 0.), axis=-1, keepdims=True)

    y_true_ = tf.boolean_mask(y_true, mask)
    y_pred_ = tf.boolean_mask(y_pred, mask)

    return keras.backend.sqrt(keras.backend.mean(keras.backend.square(y_pred_ - y_true_)))

hl1 = 125
hl2 = 125
window_len = 30
n_features = 50
batch_size = 128

optimizer = keras.optimizers.Adam(lr=0.0001, beta_1=0.9, beta_2=0.999, epsilon=None, decay=0., amsgrad=False)
dropout = 0.
input_ = keras.layers.Input(
        shape=(window_len, n_features)
    )

# masking is to make sure the model doesn't fit the zero paddings
masking = keras.layers.Masking(mask_value=0.0)(input_)

# hidden layer 1 with he_normal initializer. 

lstm_h1 = keras.layers.LSTM(hl1, dropout=dropout, kernel_initializer='he_normal',
                        return_sequences=True)(masking)

# hidden layer 2
lstm_h2 = keras.layers.LSTM(hl2, dropout=dropout, kernel_initializer='he_normal',
                        return_sequences=True)(lstm_h1)

 # dense output layer of single output
 out1 = keras.layers.Dense(1, activation='linear',name='out1')(lstm_h2)

 out2 = keras.layers.Dense(1, activation='linear', name='out2')(lstm_h2)

 out3 = keras.layers.Dense(1, activation='linear', name='out3')(lstm_h2)

 model = keras.models.Model(inputs=input_, outputs=[out1, out2, out3])


 pi = [rmse]

 n_gpus = len(get_available_gpus())

 if n_gpus > 1:
    print("Using Multiple GPU's ...")
    parallel_model = multi_gpu_model(model, gpus=n_gpus)

else:
    print("Using Single GPU ...")
    parallel_model = model


 parallel_model.compile(loss=masked_mse, optimizer=optimizer, metrics=pi)
 parallel_model.summary()

 checkpoint = keras.callbacks.ModelCheckpoint(
             file_name+".hdf5", monitor='val_loss', verbose=1, 
             save_best_only=True, mode='min', period=10,
             save_weights_only=True)

 save_history = keras.callbacks.CSVLogger(file_name+".csv", append=True)

 callbacks_list = [checkpoint, save_history]

 y_train_reshaped = list(reshape_test(window_len, y_train))
 parallel_model.fit(
            x_train,

            {
                 'out1': y_train_reshaped[0],
                 'out2': y_train_reshaped[1],
                 'out3': y_train_reshaped[2],
            },

            epochs=epochs,
            batch_size=batch_size,
            verbose=0,
            shuffle='batch',
            validation_data=(x_test, list(reshape_test(window_len,y_test))),
            callbacks=callbacks_list,

        )