shape=hermitian can be used if the matrices are going to be complex. This forces real eigenvalues for a floating point calculation.
You are worried that two zero eigenvalues not being exactly the same, and suggest that the one that is not exactly zero has the wrong symmetry for its eigenvector. But from a numerical point of view at Digits=15, these are both zero. I think the eigenvector algorithm is robust, and the symmetry issue is not related to the numerical roundoff in the eigenvalue, if that is what you are suggesting. You can take the Matrix with integer (or exact rational) entries and use NullSpace(M) to see an exact pair of eigenvectors for the zero eigenvalues. Your <0,-1,0,-1> and <1,0,-1,0> are an exact pair.
Probably you know this already, but the non-degenerate eigenvalues have unique eigenvectors up to a scalar multiplication, so the eigenvectors for E and -E will automatically have the required symmetry if E and -E are non-degenerate. Then you can just sort the eigenvectors to have them next to each other (and scale them) if that is what you want.
Once you have degeneracy, then mathematically there is no preferred form for the eigenvectors. So what you want to do is find a linear combination of the ones you get in some prescribed form. I don't see a simple way to do this. On the other hand, degeneracy usually arises from some sort of symmetry, so knowing more about the matrix structure might help you, e.g., if you can partition them into related blocks.
I'm still confused over whether you think the zero eigenvectors can alsways have the <u,v> <v,u> symmetry as I think you originally proposed, or now you think <u,u>, <v,v>?