I don’t have to use optimizer.zero_grad() before loss.backward(), I can also use it after loss.backward(), but I think it works the same no matter which version I used, each time when the train function is called, so for each batch the gradient is set to zero.
I made some changes, now I accumulate decoder outputs, and then calculate loss on entire batch.
 for input_seq, target_seq, input_lengths, target_lengths in train_loader:
    # Run a training iteration with batch

    output = train(input_seq, target_seq, input_lengths, target_lengths, encoder, decoder, embedding,
      encoder_optimizer, decoder_optimizer, clip, teacher_forcing_ratio, criterion)
    
    #print(output.contiguous().view(-1, len(voc.word_count)).shape)
    #print(target_seq.contiguous().view(-1).shape)

    loss = criterion(output.view(-1, len(voc.word_count)), 
                     target_seq.contiguous().view(-1))
    loss = loss / target_seq.size(0)

    loss.backward()          
    encoder_optimizer.step()
    decoder_optimizer.step()
    encoder_optimizer.zero_grad()
    decoder_optimizer.zero_grad()

`
In train function loop through words in given sequences:
 #Forward propagate through decoder using SOS tokens
decoder_output, decoder_hidden = decoder(decoder_input, decoder_hidden, encoder_outputs)

output_vectors.append(decoder_output)

#Forward batch of sequences through decoder one time step at a time
if use_teacher_forcing:
    for t in range(seq_len-1):
        
        # Teacher forcing: current target word is the next input
        decoder_input = target_seq[:, t].view(-1,1) # (batch_size, 1)
        
        # Forward propagate through decoder
        decoder_output, decoder_hidden = decoder(decoder_input, decoder_hidden, encoder_outputs)
        
        output_vectors.append(decoder_output)
 return torch.cat(output_vectors, dim=1)

`
When calculation loss decoders inputs correspond to targets in following manner:
decoder_input (train on) SOS->1->2->3->4
targets                                   1->2->3->4->END
In loss function I used following shapes:
loss = criterion(output.view(-1, len(voc.word_count)),
target_seq.contiguous().view(-1))
output shape: [batch_size * seq_len, vocab_size]
target shape: [batch_size * seq_len]
It is totally different approach, but loss still fluctuates at the same level.
Any ideas what can be wrong here?
Best regards.
Radek