Singular Values of Matrix as Optimization Problem Assume that $A$ is a positive semidefinite symmetric matrix. 
It is known that $$\max_{||y||\leq1} \quad y^TAy$$
Has an analytical solution which is the maximum eigenvalue of $A$. 
This isn't hard to derive either. We can do so by writing it as a constrained minimization problem and solving using KKT. 
\begin{align}
\min_{||y||\leq1} &\quad -y^TAy\\
L(y,\lambda)&=-y^TAy-\lambda(1-y^Ty)\\
\nabla_yL(y,\lambda)&=-Ay+\lambda y\\
\mbox{Setting }\nabla_yL(y,\lambda)&=0\\
Ay&=\lambda y
\end{align}
It is easy to see that the $y$ which satisfy the KKT conditions are the eigenvectors of $A$ and the highest eigenvalue gives the best objective function value. 

My question is, if now, $A$ is a rectangular matrix, is the largest singular value of $A$, the solution to the problem $$\max_{\|x\|\leq1,\|y\|\leq1} \quad x^TAy$$
where $A\in \Re^{m\times n}, x\in \Re^{m\times 1}, y\in \Re^{n\times 1}$
 A: The singular values of $A$ are the eigenvalues of the symmetric matrix
$$\begin{bmatrix}0&A\\A^T&0\end{bmatrix}$$
Following your first insight, this would mean that the largest eigenvalue can be characterized by
$$\max_{\|x\|^2+\|y\|^2\le 1} 2x^TAy$$
The set of vectors $(x,y)$ with  $\|x\|\le 1$ and $\|y\|\le1$ is strictly larger than the unit ball of the combined dimensions, so that your second formula will give the wrong value.

Let's check what the torus variant comes up to. The Lagrangian is 
$$L(x,y,λ,μ)=x^TAy+\fracλ2(1-x^Tx)+\fracμ2(1-y^Ty)$$
resulting in the KKT equations $Ay=μx$, $A^Tx=λy$, and the original norm conditions, so that 
$$μ=μx^Tx=x^TAy=λy^Ty=λ\text{ and } AA^Tx=λμx=λ^2x\text{ or } A^TAy=λμy=λ^2y.$$
Thus also resulting in the correct singular value and singular vectors.
A: Yes, that's true. It follows easily from two simple facts:
1) For $\|x\|=\|y\|=1$, $x^TAy$ is bounded from above by the maximal singular value of $A$.
If $\|x\|=\|y\|=1$, then $x^TAy\leq\|x\|\|Ay\|\leq\|x\|\|A\|\|y\|=\|A\|=\sigma_1$.
2) $x^TAy=\sigma_1$ for some $x$ and $y$ such that $\|x\|=\|y\|=1$.
Let $A=USV^T$ be the SVD with the singular values in the descending order on the diagonal of $S$ and let $e_1$ be the first column of the identity matrix of an appropriate dimension. Then taking $x=Ue_1$ and $y=Ve_1$ gives
$x^TAy=e_1^TU^TUSV^TVe_1=e_1^TSe_1=\sigma_1$.
