Very interesting approach. So, the loss function is the sum of these two terms? And the
   ( (x_i**2 + y_i**2) - 1 )**2,

term refers only to the prediction?
EDIT: Supposing I have a vector of 1024 outputs with (x_i, y_i) pairs I think what you describe should look like this:
def Cosine(output, target):

    '''

    Cosine loss function. Penalizes the area out of the unit circle and wraps the output around 2pi.

    '''    

    loss_1 = (np.add.reduceat(output ** 2, np.arange(0, len(output), 2)) - 1) ** 2 # penalize if output is outside the unit circle

    loss_1 = torch.from_numpy(loss_1)

    loss_2 = 1. - torch.cos(output - target) # wrap loss around 2pi

    return torch.mean(loss_1) + torch.mean(loss_2)

loss_1 calculation just takes the squares of all elements of the output array and sums them up in pairs pair two, then subtract by 1 and take the square.