Okay, another update. I found the problem. It was indeed the error function… I was using the standard reduction mode which averages the losses. Using this instead solves it:
lossfct = torch.nn.L1Loss(reduction='none')
[...]
for l in loss:
    l.mean().backward(retain_graph=True)

instead of
lossfct = torch.nn.L1Loss()
[...]
loss.backward()

I am glad I finally found a solution. But I am still interested in the cause of the differences. I am assuming that when I use the default reduction, there is some error due to the limited size of floats/doubles which then gets also backpropagated. So the reason is some kind of rounding error.
However, the result between CPU and GPU still differs. What could be the reason for this?
And another question:
Can the retain_graph=True be somehow avoided in this case? It increases memory usage
EDIT: Yes, it can by using reduction=‘sum’ instead of reduction=‘none’ or reduction=‘mean’
The mean operation is causing the differences