There are two ways:

You can parameterize your weight tensor, e.g. eigen decomposition, or just simply copying the entries reflected by the diagonal.
You can somehow project the weight matrix to your constrained subspace where it is symmetrical after each optimization step. However, if the gradients are not symmetrical (i.e. not satisfying the property you said above), then projection may or may not work well as this step is hidden away from the optimizing algorithm.

It seems that (1) is better, but you need to try to know if it works well in practical. By the way, by using first order optimization algorithms as we mostly do in deep learning, unless your model is fairly simple, you are not guaranteed to achieve lowest loss. Often even local minima is not gauranteed.