Setting:
Given a set of $n\times n$ matrices $A_i$, I would like to find a linear combination of these matrices $Q = \sum_i A_i x_i$ with $x_i$ a set of complex numbers, such that $Q$ is unitary: $Q^{\dagger} Q = 1$. This problem is equivalent to solving a system of quadratic equations over real numbers. As far as I understand, there is no general and efficient way to solve such systems of equations, and black-box algorithms, such as Gröbner basis, struggle with systems of around 10 variables.
Question:
Does the particular structure of this system of equations make it easier than a generic one, and can this be utilized in order to speed up the calculation.
Motivation:
This problem arises in many contexts. For example, it is related to search of symmetry of a quantum Hamiltonian. Hamiltonian $H$, a finite Hermitian matrix has a symmetry if it commutes with some unitary matrix $U$ other than identity. It is very easy to construct the linear space to which $U$ should belong, it is given by the kernel of the system of linear equations $H A - A H = 0$. However the next step requires verifying whether this linear space contains unitary matrices other than identity.
Secondary question:
Given that it seems unlikely that there is an easy answer to the main question, I would also like to ask whether there are known classes of systems of quadratic equations that are quickly solvable numerically.