Tensorflow: Can't overfit training data with b

I coded a small RNN network with Tensorflow to return the total energy consumption given some parameters. There seem to be a problem in my code. It can't overfit the training data when I use a batch size > 1 (even with only 4 samples!). In the code below, the loss value reaches 0 when I set BatchSize to 1. However, by setting BatchSize to 2, the network fails to overfit and the loss value goes toward 12.500000 and gets stuck there forever.

I suspect this has something to do with LSTM states. I get the same problem if I don't update the state with each iteration. Or maybe the cost function? A help is appreciated. Thanks.

import tensorflow as tf
import numpy as np
import os

from utils import loadData

Epochs = 10000
LearningRate = 0.0001
MaxGradNorm = 5

SeqLen = 1
NChannels = 28
NClasses = 1

NLayers = 2
NUnits = 256

BatchSize = 1
NumSamples = 4
#################################################################

trainingFile = "./training.dat"

X_values, Y_values = loadData(trainingFile, SeqLen, NumSamples)

X = tf.placeholder(tf.float32, [BatchSize, SeqLen, NChannels], name='inputs')

Y = tf.placeholder(tf.float32, [BatchSize, SeqLen, NClasses], name='labels')

keep_prob = tf.placeholder(tf.float32, name='keep')

initializer = tf.contrib.layers.xavier_initializer()

Xin = tf.unstack(tf.transpose(X, perm=[1, 0, 2]))

lstm_layers = []

for i in range(NLayers):

    lstm_layer = tf.nn.rnn_cell.LSTMCell(num_units=NUnits, initializer=initializer, use_peepholes=True, state_is_tuple=True)

    dropout_layer = tf.contrib.rnn.DropoutWrapper(lstm_layer, output_keep_prob=keep_prob)

    #[LSTM ---> DROPOUT] ---> [LSTM ---> DROPOUT] ---> etc...
    lstm_layers.append(dropout_layer)   

rnn = tf.nn.rnn_cell.MultiRNNCell(lstm_layers, state_is_tuple=True)

initial_state = rnn.zero_state(BatchSize, tf.float32)

outputs, final_state = tf.nn.static_rnn(rnn, Xin, dtype=tf.float32, initial_state=initial_state)

outputs = tf.transpose(outputs, [1,0,2])
outputs = tf.reshape(outputs, [-1, NUnits])

weight = tf.Variable(tf.truncated_normal([NUnits, NClasses]))
bias = tf.Variable(tf.constant(0.1, shape=[NClasses]))
prediction = tf.matmul(outputs, weight) + bias
prediction = tf.reshape(prediction, [BatchSize, SeqLen, NClasses])

cost = tf.reduce_sum(tf.pow(tf.subtract(prediction, Y), 2)) / (2 * BatchSize)

tvars = tf.trainable_variables()

grad, _ = tf.clip_by_global_norm(tf.gradients(cost, tvars), MaxGradNorm)

optimizer = tf.train.AdamOptimizer(learning_rate = LearningRate)

train_step = optimizer.apply_gradients(zip(grad, tvars))

sess = tf.Session()

sess.run(tf.global_variables_initializer())

iteration = 1

for e in range(0, Epochs):

    train_loss = []

    state = sess.run(initial_state)

    for i in xrange(0, len(X_values), BatchSize):
        x = X_values[i:i + BatchSize]
        y = Y_values[i:i + BatchSize]

    y = np.expand_dims(y, 2)

        feed = {X : x, Y : y, keep_prob : 1.0, initial_state : state}

        _ , loss, state, pred = sess.run([train_step, cost, final_state, prediction], feed_dict = feed)

        train_loss.append(loss)

        iteration += 1

    print("Epoch: {}/{}".format(e, Epochs), "Iteration: {:d}".format(iteration), "Train average rmse: {:6f}".format(np.mean(train_loss)))