This isn't pretty, but it's direct. We proceed by the usual method of finding the critical points of a function followed by the critical points on the boundary. The maximum occurs at once of these points.
To find the critical points of the objective function, differentiate:
$$ df(x) = (1/2) (dx)^T A x + (1/2) x^T A dx + b^T dx = (dx)^T (Ax + b) $$
so the critical points are solutions to $Ax = -b$.
Using Lagrange multipliers to find the critical points on the boundary of your constrant $x^T x = 1$ says they occur at solutions $(x, \mu)$ to the system of equations
$$ (A + 2 \mu I) x = -b \qquad \qquad \qquad x^T x = 1$$
Most solutions can be obtained by solving the equation
$$ || (A + 2 \mu I)^{-1} b ||_2^2 = 1 $$
for $\mu$. However, we must consider separately the points $\mu$ for which $-2\mu$ is an eigenvalue of $A$.
For each of these $\mu$, obtain the general solution to
$$ (A + 2 \mu I) x = -b $$
and and search the solution space for unit vectors. Typically there won't be any solutions at all! But if there are, you have to search the solution space for unit vectors. This is just solving a polynomial equation in $n$ variables, where $n$ is the multiplicty of the eigenvalue $-2 \mu$.
I get essentially the same thing from the eigen-decomposition. Write everything in terms of an orthonormal basis of eigenvectors. Then, you are trying to maximize
$$ (1/2) \sum_i x_i^2 \lambda_i + x_i b_i$$
subject to the constraint
$$ \sum_i x_i^2 = 1 $$
The criticial points of the objective function are the solutions to
$$ x_i \lambda_i + b_i = 0$$
and Lagrange multipliers gives us a system of equations
$$ \sum_i x_i^2 = 1 \qquad \qquad
x_i( \lambda_i + 2 \mu) + b_i = 0 $$
and the solution method to this is effectively the same as the matrix version.
Let $\lambda_{\min}$ be the minimum eigenvalue of $Q$.
The dual function is
\begin{align}
g(u) &= \inf_x L(x,u) \\
&= \begin{cases}
-\infty \quad \text{if } u < -\lambda_{\min} , \\
-\frac12 u \quad \text{otherwise}.
\end{cases}
\end{align}
The dual problem is
\begin{align}
\operatorname*{maximize}_{u} & \quad -\frac12 u \\
\text{subject to} & \quad u \geq -\lambda_{\min}.
\end{align}
The optimal value for the dual problem is
\begin{equation*}
d^\star = \frac12 \lambda_{\min}.
\end{equation*}
And this is equal to the primal optimal value, according to the standard variational characterization of the eigenvalues of $Q$.
Note: To derive the dual function, I used the following facts.
- The minimum eigenvalue of $Q + u I$ is $\lambda_\min + u$.
- Let $M \in \mathbb R^{n \times n}$ be symmetric. The quadratic function $h(x) = \frac12 x^T M x$ is unbounded below if the minimum eigenvalue of $M$ is negative. Otherwise, the minimum value of $M$ is $0$.
Best Answer
This is really more of a comment than a complete answer, but in case you are not aware, a slightly simpler version of your quadratic program can be reformulated as,
\begin{align} \text{minimize}\;\;&\langle A,X\rangle\\[4pt] \text{subject to}\;\;&\langle B,X\rangle=0\\[1pt] &\text{rank}(X)=1\\[1pt] &X\succeq 0 \end{align}
where the final two constraints ensure that $x$ may be recovered from $X$ via its outer product (i.e. $X=xx^T$).
Your actual problem can be handled in a similar fashion by first using a process called homogenization.
The non-convex aspect has now been isolated in the rank constraint, which, if removed, will yield a semidefinite relaxation of your quadratic program, this relaxation is canonical, and people have studied it in detail (I have not however).