I came up with a solution that avoids hard constraints by increasing the loss if the weights are not symetrical.
def sym_loss(w, alpha=0.1, dims=(0,)):
    nb_w = w.shape[0]*w.shape[1]
    return alpha * torch.pow(w - w.flip(dims=dims), 2).sum().sqrt()/nb_w

Here, w corresponds to the weights of the convolutional layer of interest. The default dims enforces an horizontal symetry. The coefficient alpha indicates the importance of this loss in the overall loss.
If you have any feedback on this approach, please comment.
Personnaly I see one draw back which is that you still have to search and update weights on both sides of the symetrical axis while you could constrain your search space more efficiently and precisely is a string constraint defined in advance. Here you won’t have a guarentee of a symetry.