Thank you for your answer.
I tried setting the Softmax to a LogSoftmax and the KL divergence error is indeed positive now. The target values look something like this
tensor([[0.000 0.000 0.149 0.000 0.000 0.000 0.851 0.000 0.000 0.000] ... (and then 2 more)]])

Which should be correct since they add up to 1, right? After chancing the Softmax to LogSoftmax however, the output values look something like this
tensor([[-2.129 -2.201 -2.1317 ....

Which is not a normal probability distribution anymore. Should anything be changed to the target values as well?
Another problem I had is in the network architecture, in the paper they use
('conv1', ConvBNReLu(3, 64, 3, padding=1)),('conv2', ConvBNReLu(64, 64, 3, padding=1)) ...

which gave me the error:
RuntimeError: Given groups=1, weight of size 32 32 3 3, expected input[3, 64, 32, 32] to have 32 channels, but got 64 channels instead


Which was solved by halving all out_channel values, however the authors of the paper do use a higher number of filters, looking at Table 7 in the paper. What am I doing wrong there?