That follows is a detailed answer to a comment of Federico Poloni.

Since the OP does not seem interested in a successive approximation resolution, I consider the problem from an algebraic point of view. In the sequel, we assume that $n\geq 2$.

Let $A=[a_{i,j}],B=[b_{i,j}]$ be generic $n\times n$ real matrices (the $(a_{i,j}),(b_{i,j})$ are parameters that are mutually transcendental over $\mathbb{Q}$) . We consider the quotient field $K=\mathbb{Q}((a_{i,j}),(b_{i,j}))$.

We consider the problem $(\mathcal{P})$: search the maximum of the function

$f:(x,y)\in \mathbb{R}^n\times \mathbb{R}^n\rightarrow (Ax+By)^T(Ax+By)$ under the conditions $(1)$ $x^Tx=y^Ty=1$.

$\textbf{Lemma 1}$. The Lagrange condition for a local extremum of $f$ under $(1)$ is

$(2)$ $A^T(Ax+By)-ux=0,B^T(Ax+By)-vy=0$ where $u,v\in\mathbb{R}$.

$\textbf{Proof}$. There are real $u,v$ s.t. for every vectors $h,k$, one has

$h^TA^T(Ax+By)+k^TB^T(Ax+By)-uh^Tx-vk^Ty=0$,

That implies $(2)$. $\square$

$\textbf{Lemma 2}$. The absolute minimum of $f$ is obtained on $U=\{(x,y);Ax+By=0,x^Tx=y^Ty=1\}$, an algebraic set of dimension $n-2$.

$\textbf{Proof}$. If $(x,y)\in U$, thet $A(x)=B(-y)\in V=A(S^{n-1})\cap B(S^{n-1})\subset \mathbb{R}^n$. When $A,B$ are in general position, $V$ is an algebraic set of dimension $2(n-1)-n=n-2$ which is diffeomorphic to $U$. $\square$

Therefore the system $\mathcal{S}=\{(1),(2),Ax+By\not= 0\}$ provides all the candidates for $\max(f)$. We go to see that the associated ideal has a Hilbert dimension equal to $0$.

$\textbf{Proposition}$. When $n\leq 6$, $\mathcal{S}$ has generically at most $2n(n+1)$ real solutions; moreover, the complexity of solving $\mathcal{S}$ is the same as the complexity of solving a polynomial of degree $n(n+1)$
with Galois group $S_{n(n+1)}$ over $K$. In particular, the generic problem $\mathcal{P}$ is non-solvable by radicals.

$\textbf{Proof}$. We use Grobner basis theory, the unknowns being $x=[x_i],y=[y_i]$; that's why I guess $n\leq 6$ (the time of calculation for $n=6$ is $2$ minutes).

Generically, the system reduces to a system in the form

$\{P_{2n(n+1)}(x_1)=0,x_i=Q_i(x_1),y_j=R_j(x_1)\}$ where $P_{2n(n+1)}$ is an even polynomial of degree $2n(n+1)$ and $Q_i,R_j$ are polynomials of degree $<2n(n+1)$, all being given explicitly by the software (note that, for $n=6$, the size of the coefficients is huge!). The parity of $P(x_1)=\tilde{P}(x_1^2)$ comes from $f(x,y)=f(-x,-y)$.

Clearly, the complexity lies entirely in the search for roots of $\tilde{P}$. To show that the Galois group of $\tilde{P}$ is $S_{n(n+1)}$, it suffices to use the "specialization theorem" that says that if we choose explicit values for the $(a_{i,j}),(b_{i,j})$, then the Galois group of the obtained $\tilde{P}_0$ is a subgroup of the Galois group associated to the generic $\tilde{P}$.

That can be done (with random choices in $\mathbb{Z}$) for $n\leq 6$. $\square$

$\textbf{Conjecture}$. The result of the above Proposition is true for every $n$.

$\textbf{Remarks}.$ 1. When we know the $O(n^2)$ candidates $(x^i,y^i)$, to obtain the required maximum, it suffices to test the associated values of $f(x^i,y^i)$, that has a total complexity in $O(n^4)$.

However, for every $n\geq 2$, the generic polynomial $\tilde{P}$ is non-solvable (by radicals). That implies that, if we randomly choose $A,B$ (the $(a_{i,j}),(b_{i,j})$ are independent and follow a normal law), then the problem $\mathcal{P}$ is non-solvable (by radicals) with probability $1$. We can calculate an approximation of the roots of $\tilde{P}$ with complexity $O(n^3)$.

Of course, there are couples $(A,B)$ s.t. $\mathcal{P}$ is solvable. For example, when $A,B\in O(n)$, $\max(f)=2^2=4$.

Convex Optimization). – Suvrit Oct 22 at 2:48