According to the source code:

  
      github.com
  
  
    tensorflow/tensorflow/blob/v2.1.0/tensorflow/python/keras/layers/recurrent_v2.py#L890-L1170 1

class LSTM(recurrent.DropoutRNNCellMixin, recurrent.LSTM):
  """Long Short-Term Memory layer - Hochreiter 1997.


  See [the Keras RNN API guide](https://www.tensorflow.org/guide/keras/rnn)
  for details about the usage of RNN API.


  Based on available runtime hardware and constraints, this layer
  will choose different implementations (cuDNN-based or pure-TensorFlow)
  to maximize the performance. If a GPU is available and all
  the arguments to the layer meet the requirement of the CuDNN kernel
  (see below for details), the layer will use a fast cuDNN implementation.


  The requirements to use the cuDNN implementation are:


  1. `activation` == `tanh`
  2. `recurrent_activation` == `sigmoid`
  3. `recurrent_dropout` == 0
  4. `unroll` is `False`
  5. `use_bias` is `True`
  6. Inputs are not masked or strictly right padded.


  This file has been truncated. show original

  
  
    
    
  
  


    initial_state: List of initial state tensors to be passed to the first
      call of the cell (optional, defaults to `None` which causes creation
      of zero-filled initial state tensors).

So the hidden state is always zero-filled just like PyTorch. Am I possibly doing something wrong during the actual training phase? I have another model in PyTorch which I train the same way, so I don’t think this is the issue, but it may be the case here.
# function: train                                                                                                                                                                                                   
#                                                                                                                                                                                                                   
# arguments: tr_gen - data loader for training data                                                                                                                                                                 
#            cv_gen - data loader for cross_val data                                                                                                                                                                
#            model - the pytorch model                                                                                                                                                                              
#            opt - the optimizer for the model                                                                                                                                                                      
#            loss_fx - the loss function used to backprop                                                                                                                                                           
#            device - cpu/gpu device used for operations                                                                                                                                                            
#                                                                                                                                                                                                                   
# return: history - acc/loss dict for training and cv                                                                                                                                                               
#                                                                                                                                                                                                                   
# This method trains on the dataset for one epoch                                                                                                                                                                   
#
def train(tr_gen, cv_gen, model, opt, loss_fx, device):

    # set data loaders to a dictionary                                                                                                                                                                              
    #                                                                                                                                                                                                               
    data_loaders = {TRAIN: tr_gen, EVAL: cv_gen}

    # set losses for each phase to a dictionary                                                                                                                                                                     
    #                                                                                                                                                                                                               
    phase_losses = {TRAIN: 0.0, EVAL: 0.0}

    # set accuracies for each phase to a dictionary                                                                                                                                                                 
    #                                                                                                                                                                                                               
    phase_acc = {TRAIN: 0.0, EVAL: 0.0}

    # time an epoch                                                                                                                                                                                                 
    #                                                                                                                                                                                                               
    start_time = time.time()

    # for each phase (train/cv)                                                                                                                                                                                     
    #                                                                                                                                                                                                               
    for phase in [TRAIN, EVAL]:

        # set the approriate mode for the model                                                                                                                                                                     
        #                                                                                                                                                                                                           
        if phase == TRAIN:
            model.train(True)
        else:
            model.eval()
        # collect the loss average over batches                                                                                                                                                                     
        #                                                                                                                                                                                                           
        running_loss = 0.0

        # collect the number of correct classifications                                                                                                                                                             
        #                                                                                                                                                                                                           
        correct = 0

        # total number of inputs to the model                                                                                                                                                                       
        #                                                                                                                                                                                                           
        total = 0

        # collect number of batches                                                                                                                                                                                 
        #                                                                                                                                                                                                           
        num_batches = len(data_loaders[phase])

        # get the batch size                                                                                                                                                                                        
        #                                                                                                                                                                                                           
        batch_size = data_loaders[phase].dataset.frames_per_minibatch

        # iterate over the dataset                                                                                                                                                                                  
        #                                                                                                                                                                                                           
        for index, batch in enumerate(data_loaders[phase]):

            # zero out the optimizer gradient                                                                                                                                                                       
            #                                          
            opt.zero_grad()

            # collect the data and labels                                                                                                                                                                           
            #                                                                                                                                                                                                       
            data = batch[0].squeeze().to(device)
            labels = batch[1].squeeze().max(dim=1)[1].to(device)
            batch_size = data.shape[0]

            # feed the data to the model                                                                                                                                                                            
            #                                                                                                                                                                                                       
            output = model(data)

            # calculate the loss                                                                                                                                                                                    
            #                                                                                                                                                                                                       
            loss = loss_fx(output, labels)
            running_loss += loss.item()

            # increment the number of input to the model                                                                                                                                                            
            #                                                                                                                                                                                                       
            total += labels.shape[0]

            # count the times where the label index == output max index                                                                                                                                             
            #                                                                                                                                                                                                       
            correct += torch.sum(labels == torch.max(output.data, dim=1)[1].to(device))

            # if we are training                                                                                                                                                                                    
            #                                       
            if phase == TRAIN:

                # perform backprop and take a step for the optimizer                                                                                                                                                
                #                                                                                                                                                                                                   
                loss.backward()
                opt.step()

        # calculate the average loss for the dataset                                                                                                                                                                
        #                                                                                                                                                                                                           
        phase_losses[phase] = running_loss / float(num_batches)

        # calculate the accuracy for all input                                                                                                                                                                      
        #                                                                                                                                                                                                           
        phase_acc[phase] = float(correct) / float(total)

    # collect the loss and accuracy                                                                                                                                                                                 
    #                                                                                                                                                                                                               
    history = {LOSS: phase_losses, ACC: phase_acc}

    # calculate the time it took to train                                                                                                                                                                           
    #                                                                                                                                                                                                               
    total_time = int(time.time() - start_time)

    # collect loss/acc for training and cv                                                                                                                                                                          
    #                                                                                                                                                                                                               
    loss = phase_losses[TRAIN]
    acc = phase_acc[TRAIN]
    val_loss = phase_losses[EVAL]
    val_acc = phase_acc[EVAL]

    # print informational message                                                                                                                                                                                   
    #                                                                                                                                                                                                               
    print(" - %ds - loss: %.4f - acc: %.4f - val_loss: %.4f - val_acc: %.4f"
          % (total_time, loss, acc, val_loss, val_acc))

    # return the history                                                                                                                                                                                            
    #                                                                                                                                                                                                               
    return history
#                                                                                                                                                                                                                   
# end of function

It should be noted that the way I train the model is by setting a learning rate initially and manually reducing it by using this function:
# function: update_lr                                                                                                                                                                                               
#                                                                                                                                                                                                                   
# arguments: optim - the adam optimizer                                                                                                                                                                             
#            lr - the new learning rate                                                                                                                                                                             
#                                                                                                                                                                                                                   
# return: none                                                                                                                                                                                                      
#                                                                                                                                                                                                                   
# This method updates the learning rate                                                                                                                                                                             
#                                                                                                                                                                                                                   
def update_lr(optim, lr):

    # update the learning rate for each param group                                                                                                                                                                 
    #                                                                                                                                                                                                               
    for param_group in optim.param_groups:
        param_group['lr'] = lr
#                                                                                                                                                                                                                   
# end of function 

The learning rate is halved each epoch once the validation error after the current epoch is at least >= min_val_error, which is just a custom set parameter. The same technique is used for the Keras model.
Thanks again for your help!