Coupling of Gaussian Random variables. Suppose I have two mean zero multivariate Gaussian random variables $X$ and $Y$ on $\mathbb R^d$, with covariance matrices $A$ and $B$. (Assumed to have full rank).
I have many choices of joint distributions such that the pair $(X,Y)\in\mathbb R^{2d}$ is Gaussian.  I want to choose one so that $X$ and $Y$ are as close to each other as possible.
So as $X-Y$ is a Gaussian I want to minimise the trace of its covariance matrix. 
I'm pretty sure the way to do this is to choose some matrix $M$ such that $B =M^T A M$ and set $Y = M^T\cdot X$.  
Intuitively, this reduces the range of values that the pair $(X,Y)$ can take, so it must reduce the variance of $X-Y$.  For my purposes this reasoning is good enough.
I already have $A$ and $B$ in the form $A= \alpha^T\alpha$, $B=\beta^T\beta$ so I think I just need to choose a unitary matrix $U$ and set $M = \alpha^{-1} U \beta$.
I've managed to forget everything I know about linear algebra. So I'm having trouble minimising over $U$, the naïve way of doing it (differentiate everything in sight and set to $0$) just produces a mess.
Is there a nice way of doing this?
 A: The trace of the covariance matrix of $X-Y$ is
$$\def\Tr{\operatorname{Tr}}\Tr\langle(X-Y)^\top(X-Y)\rangle=\langle\Tr X^\top X+\Tr Y^\top Y-2\Tr X^\top Y\rangle\;,$$
where angled brackets denote the expected value. The only part that varies with $U$ is $\langle\Tr X^\top Y\rangle$. 
Since Michael has now confirmed that there's an error in your ansatz, I'll write out the calculation again in what I think is a correct form.
In order to leave most of your setup of $\alpha$, $\beta$, $A$, $B$ and $U$ intact, I'll just replace $X=MY$ by $M^\top X=Y$. Then we have
\begin{align}
\langle\Tr X^\top Y\rangle
&=
\langle\Tr X^\top M^\top X\rangle
\\
&=
\langle\Tr X^\top\beta^\top U^\top\alpha^{-1\top} X\rangle
\\
&=
\Tr\alpha^{-1\top}\langle XX^\top\rangle\beta^\top U^\top
\\
&=
\Tr\alpha^{-1\top}A\beta^\top U^\top
\\
&=
\Tr\alpha^{-1\top}\alpha^\top\alpha\beta^\top U^\top
\\
&=
\Tr\alpha\beta^\top U^\top
\;.
\end{align}
Thus the optimal $U$ is the unitary matrix closest to $\alpha\beta^\top$, which you can compute using singular value decomposition. As clarified in an exchange of comments between Michael and myself, the minimization of $\lVert R-A\rVert_F$ mentioned in that Wikipedia article is equivalent to the maximization of $\operatorname{Tr}A^\top R$, since $\lVert R-A\rVert_F^2=\operatorname{Tr}(R^\top R-R^\top A-A^\top R+A^\top A)$, where the first and last terms are constant and the other two both yield $-\operatorname{Tr}A^\top R$.
