If we didn’t worry about DCP compliance, the constraints could be written
U(i,:)*A*U(i,:)' == 0
However, those constraints are all non-convex as written. Because
U is binary,
U(i,:)*A*U(i,:)' can be linearized. To do so, you can do the following, which should work, but I’m not sure is the best or fastest way of doing this.
Introduce another binary variable for every product term
U(i,j)*(U(k,l). Call that
B(i,j,k,l).for each instance of
U(i,j)*U(k,l). Then add the constraints
B(i,j,k,l) <= U(i,j)
B(i,j,k,l) <= U(k,l)
B(i,j,k,l) >= U(i,j) + U(k,l) - 1
I think you really only need
(M*N)*(M*N -1)/2 of these variables due to symmetry and that square terms (i.e., of the form
U*i,j)*U(i,j) ) can be replaced by the non-squared variable without adding variables or constraints.
Anyhow, it might get a bit messy, but I leave you to deal with the mess and to determine any efficient vectorized ways of implementing this. Converting from a Binary QP to a MILP is a standard thing, so that’s basically what you have to do (pretend the left hand side of the quadratic equality constraint is the objective function). Keep in mind though that your optimization variable is a matrix, not a vector.
Edit: I leave you to determine whether there is a better way of dealing with the requirement stated at the end of the edited version of your most recent post. I wouldn’t be surprised if there were a much nicer way, requiring many fewer variables and constraints.