If you can translate that loss function into a python expression using torch compatible math ops, then you can do
loss = that_expression
# and then you can get the gradients and update as usual
loss.backward()
optimizer.step()

No need to implement the backward pass explicitly.
Some torch compatible math ops include…
torch.sin(var_theta)
var_a * var_b
var_a + var_b
var_a - var_b
torch.sum(variable, dim=dimension_to_sum_over)

I do not think that having an infinite number of possible solutions is a problem. I would not constrain theta at all, just let it find the nearest solution. That is, unless you have a good reason for wanting to constrain the search to only small rotations of the input.