I doubt the problem arises from loss computation;
if I change the linear layer input to output[:, -1, :] (only use the last time-step output),
and compute loss as loss = loss_function(log_prob, y[:, -1]) (only use the last time-step label),
I can get the expected same result from the two inference methods.
But I would like to compute accumulated loss of every time-step instead of just the last one in a sequence.
Though feeding in a single time-step and computing the loss iteratively throughout a sequence (as follows) works, it is very time consuming during training (about four times slower as opposed to feeding in a sequence as a whole), which is why I consider the training method as shown in my post.
hidden = nodel.init_hidden()
loss = 0
for t in range(seq_len):
    log_prob, hidden = model(x_onehot[:, t, :], hidden)
    loss_t = loss_function(log_prob, y[:, t])
    loss += loss_t
loss.backward()

Could someone help me please? On efficient training of LSTMs which accumulates loss from every time-step, and behaves as expected during inference; or perhaps I am missing something.