I think there is an issue with the solution proposed here: we are fixing the model to accept only a batch size of 20, and we are learning separate values of h0 for each of the 20 batch elements.
This is how I solved it:
def __init__(...):
    # learn initial hidden state (h0, c0)
    # NOTE: use "2*" only if bidirectional
    h0 = torch.zeros(2*n_layers, 1, self.hidden_dim).to(device)
    c0 = torch.zeros(2*n_layers, 1, self.hidden_dim).to(device)
    nn.init.xavier_normal_(h0, gain=nn.init.calculate_gain('relu'))
    nn.init.xavier_normal_(c0, gain=nn.init.calculate_gain('relu'))
    self.h0 = nn.Parameter(h0, requires_grad=True)  # Parameter() to update weights
    self.c0 = nn.Parameter(c0, requires_grad=True)
    # initialize the other stuff...

def forward(self, seqs, lengths):
    packed = nn.utils.rnn.pack_padded_sequence(seqs, lengths)
    batch_size = seqs.shape[1]
    # repeat hid0 for batch while keeping the gradient update working for whole batch
    outputs, (hidden, cell) = self.lstm(packed, (self.h0.repeat(1, batch_size, 1),
                                                 self.c0.repeat(1, batch_size, 1)))
    # do more stuff / return ...

I also tried using expand() instead of repeat() to avoid creating copies. However these tensors need to be contiguous, so I think we can’t avoid the memory allocation (note that expand().contiguous() also creates a new copy). Both approaches seem to compute the gradients correctly, by considering the whole batch.
I’m a bit surprised that there isn’t a simpler way to do this (e.g. by providing the initial state h0 to the LSTM separately, for batch_size=1, and then feeding the whole batch without providing h0). If someone knows a cleaner solution please let me know.