I'm interested in
$$\sup_{x, y} \left\| A x + B y\right\|_2$$ subject to $$\left\|x\right\|_2 = \left\|y\right\|_2 = 1$$
where $A$, $B$ and $x$, $y$ are real matrices and vectors, respectively, of compatible sizes. (Generalizations to more than two terms are also interesting)
Does this problem have a name in the literature? (If not, does it reduce to a well-known one?)
Note that this is a stronger constraint than $\left\|\begin{bmatrix}x\\y\end{bmatrix}\right\|_2 = \sqrt 2$
The sets $\{Ax : \|x\|=1\}$ and $\{By : \|y\|=1\}$ are ellipsoids. Hence the set $\{Ax+By : \|x\|=\|y\|=1\}$ is the Minkowski sum of two ellipsoids. Googling for these terms returned this paper which may be interesting to you (although browsing through it I did not find an immediate solution to your problem).
Yan, Yan; Chirikjian, Gregory S., Closed-form characterization of the Minkowski sum and difference of two ellipsoids, Geom. Dedicata 177, 103-128 (2015). ZBL1321.65033.
I do not know if there is a name for this problem. However, we can view this in terms of the usual operator norm.
Given $A \in \mathrm{Mat}_{k \times n}(\mathbb{R}), B \in \mathrm{Mat}_{k \times m}(\mathbb{R})$, consider for $x \in \mathbb{R}^n, y \in \mathbb{R}^m$:
\begin{align*} \underset{\substack{ \lvert \lvert x \rvert \rvert = 1 \\ \lvert \lvert y \rvert \rvert = 1}}{\mathrm{sup}} \lvert \lvert Ax + By \rvert \rvert &= \underset{\substack{ \lvert \lvert x \rvert \rvert \leq 1 \\ \lvert \lvert y \rvert \rvert \leq 1}}{\mathrm{sup}} \lvert \lvert Ax + By \rvert \rvert \\ &= \underset{\substack{ \lvert \lvert x \rvert \rvert \leq 1 \\ \lvert \lvert y \rvert \rvert \leq 1}}{\mathrm{sup}} \, \left \lvert \left \lvert \begin{pmatrix}A & B\end{pmatrix} \begin{pmatrix}x \\ y\end{pmatrix} \right \rvert \right\rvert \\ &= \underset{\left \lvert \left \lvert \begin{pmatrix} x \\ y \end{pmatrix} \right \rvert \right \rvert_\infty = 1}{\mathrm{sup}} \left \lvert \left \lvert \begin{pmatrix}A & B\end{pmatrix} \begin{pmatrix}x \\ y\end{pmatrix} \right \rvert \right\rvert, \end{align*} where $\left \lvert \left \lvert \begin{pmatrix} x \\ y \end{pmatrix} \right \rvert \right \rvert_\infty := \max(\lvert \lvert x \rvert \rvert, \lvert \lvert y \rvert \rvert)$ is the supremum norm from viewing $\mathbb{R}^{n+m}$ as the product $\mathbb{R}^n \times \mathbb{R}^m$. Observe then that the last expression is $\lvert \lvert (\begin{smallmatrix}A & B\end{smallmatrix}) \rvert \rvert_{\mathrm{op}}$, where $\lvert \lvert \cdot \rvert \rvert_{\mathrm{op}}$ is the operator norm induced from $\lvert \lvert \cdot \rvert \rvert_\infty$ on $\mathbb{R}^{n+m}$ and $\lvert \lvert \cdot \rvert \rvert$ on $\mathbb{R}^k$.
Also note that this characterization $\underset{\lvert \lvert x_i \rvert \rvert = 1}{\mathrm{sup}} \, \lvert \lvert \sum_i A_i x_i \rvert \rvert = \lvert \lvert (\begin{smallmatrix} A_1 & \ldots & A_\ell \end{smallmatrix}) \rvert \rvert_{\mathrm{op}}$ works for any finite sum.
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$.
I write a new post (again) to show how the great idea of Suvrit (cf. above) allows to calculate the required maximum without using any algorithm; I don't calculate the $(x,y)$ that realizes the maximum.
Let $A,B\in M_n(\mathbb{R})$; we search the maximum of the function
$f:(x,y)∈\mathbb{R}^n\times\mathbb{R}^n→(Ax+By)^T(Ax+By)$ under the conditions (1) $x^Tx=y^Ty=1$.
Let $\alpha=\max(f)$; then one has the implication linking the following $3$ quadratic functions of $[x,y]^T$
$f_1=1-x^Tx\geq 0,f_2=1-y^Ty\geq 0$ $\implies$ $g=\alpha-(Ax+By)^T(Ax+By)\geq 0$.
According to the S-lemma, there exist $u,v\geq 0$ s.t., for every $x,y$, one has $g\geq uf_1+vf_2$.
Let $M_{u,v}=\begin{pmatrix}uI_n-A^TA&-A^TB\\-B^TA&vI_n-B^TB\end{pmatrix}$. The above condition is equivalent to
for every $x,y$, one has $[x^T,y^T]M_{u,v}[x,y]^T\geq u+v-\alpha$.
Necessarily, the symmetric matrix $M_{u,v}$ is $\geq 0$ and, in particular, $u\geq \rho(A^TA),v\geq \rho(B^TB)$. Moreover $x=y=0$ implies that $\alpha\geq u+v$. We search the smallest of the $\alpha$'s realizing these conditons; then $\alpha=u_0+v_0$ where $u_0+v_0$ is the minimal $u+v$ s.t. $M_{u,v}\geq 0$.
Of course, when $u,v>0$ are great enough, $M_{u,v}>0$. Then we seek $u+v$ minimal s.t. $M_{u,v}\geq 0$ and $\det(M_{u,v})=0$, that is, we search the maximal $a$ s.t. the hyperplane $u+v=a$ is tangent to the hypersurface $\det(M_{u,v})=0$. Many experiments "show" that the method works. Clearly, a rigorous proof will be welcome (in particular, in relation to the convexity of the function $\det(M_{u,v})=0$).
$\textbf{Conclusion}$. We obtain $a=\max(f)$ as follows.
i) Calculate the polynomial $g(u,a)=\det(M_{u,-u+a})$.
ii) Calculate the polynomial $h(a)=discrim(g(u,a),u)$, the discriminant of $g$ w.r.t. $u$.
iii) The required $\max$ is the greatest root of $h$.
