The purpose of this answer is to explain the exact mathematical relationship between linear discriminant analysis (LDA) and multivariate linear regression (MLR). It will turn out that the correct framework is provided by reduced rank regression (RRR).
We will show that LDA is equivalent to RRR of the whitened class indicator matrix on the data matrix.
###Notation
Let $\newcommand{\X}{\mathbf X}\X$ be the $n\times d$ matrix with data points $\newcommand{\x}{\mathbf x}\x_i$ in rows and variables in columns. Each point belongs to one of the $k$ classes, or groups. Point $\x_i$ belongs to class number $g(i)$.
Let $\newcommand{\G}{\mathbf G}\G$ be the $n \times k$ indicator matrix encoding group membership as follows: $G_{ij}=1$ if $\x_i$ belongs to class $j$, and $G_{ij}=0$ otherwise. There are $n_j$ data points in class $j$; of course $\sum n_j = n$.
We assume that the data are centered and so the global mean is equal to zero, $\newcommand{\bmu}{\boldsymbol \mu}\bmu=0$. Let $\bmu_j$ be the mean of class $j$.
LDA
The total scatter matrix $\newcommand{\C}{\mathbf C}\C=\X^\top \X$ can be decomposed into the sum of between-class and within-class scatter matrices defined as follows: \begin{align} \C_b &= \sum_j n_j \bmu_j \bmu_j^\top \\ \C_w &= \sum(\x_i - \bmu_{g(i)})(\x_i - \bmu_{g(i)})^\top. \end{align} One can verify that $\C = \C_b + \C_w$. LDA searches for discriminant axes that have maximal between-group variance and minimal within-group variance of the projection. Specifically, first discriminant axis is the unit vector $\newcommand{\w}{\mathbf w}\w$ maximizing $\w^\top \C_b \w / (\w^\top \C_w \w)$, and the first $p$ discriminant axes stacked together into a matrix $\newcommand{\W}{\mathbf W}\W$ should maximize the trace $$\DeclareMathOperator{\tr}{tr} L_\mathrm{LDA}=\tr\left(\W^\top \C_b \W (\W^\top \C_w \W)^{-1}\right).$$
Assuming that $\C_w$ is full rank, LDA solution $\W_\mathrm{LDA}$ is the matrix of eigenvectors of $\C_w^{-1} \C_b$ (ordered by the eigenvalues in the decreasing order).
This was the usual story. Now let us make two important observations.
First, within-class scatter matrix can be replaced by the total scatter matrix (ultimately because maximizing $b/w$ is equivalent to maximizing $b/(b+w)$), and indeed, it is easy to see that $\C^{-1} \C_b$ has the same eigenvectors.
Second, the between-class scatter matrix can be expressed via the group membership matrix defined above. Indeed, $\G^\top \X$ is the matrix of group sums. To get the matrix of group means, it should be multiplied by a diagonal matrix with $n_j$ on the diagonal; it's given by $\G^\top \G$. Hence, the matrix of group means is $(\G^\top \G)^{-1}\G^\top \X$ (sapienti will notice that it's a regression formula). To get $\C_b$ we need to take its scatter matrix, weighted by the same diagonal matrix, obtaining $$\C_b = \X^\top \G (\G^\top \G)^{-1}\G^\top \X.$$ If all $n_j$ are identical and equal to $m$ ("balanced dataset"), then this expression simplifies to $\X^\top \G \G^\top \X / m$.
We can define normalized indicator matrix $\newcommand{\tG}{\widetilde {\mathbf G}}\tG$ as having $1/\sqrt{n_j}$ where $\G$ has $1$. Then for both, balanced and unbalanced datasets, the expression is simply $\C_b = \X^\top \tG \tG^\top \X$. Note that $\tG$ is, up to a constant factor, the whitened indicator matrix: $\tG = \G(\G^\top \G)^{-1/2}$.
Regression
For simplicity, we will start with the case of a balanced dataset.
Consider linear regression of $\G$ on $\X$. It finds $\newcommand{\B}{\mathbf B}\B$ minimizing $\| \G - \X \B\|^2$. Reduced rank regression does the same under the constraint that $\B$ should be of the given rank $p$. If so, then $\B$ can be written as $\newcommand{\D}{\mathbf D} \newcommand{\F}{\mathbf F} \B=\D\F^\top$ with both $\D$ and $\F$ having $p$ columns. One can show that the rank two solution can be obtained from the rank solution by keeping the first column and adding an extra column, etc.
To establish the connection between LDA and linear regression, we will prove that $\D$ coincides with $\W_\mathrm{LDA}$.
The proof is straightforward. For the given $\D$, optimal $\F$ can be found via regression: $\F^\top = (\D^\top \X^\top \X \D)^{-1} \D^\top \X^\top \G$. Plugging this into the loss function, we get $$\| \G - \X \D (\D^\top \X^\top \X \D)^{-1} \D^\top \X^\top \G\|^2,$$ which can be written as trace using the identity $\|\mathbf A\|^2=\mathrm{tr}(\mathbf A \mathbf A^\top)$. After easy manipulations we get that the regression is equivalent to maximizing (!) the following scary trace: $$\tr\left(\D^\top \X^\top \G \G^\top \X \D (\D^\top \X^\top \X \D)^{-1}\right),$$ which is actually nothing else than $$\ldots = \tr\left(\D^\top \C_b \D (\D^\top \C \D)^{-1}\right)/m \sim L_\mathrm{LDA}.$$
This finishes the proof. For unbalanced datasets we need to replace $\G$ with $\tG$.
One can similarly show that adding ridge regularization to the reduced rank regression is equivalent to the regularized LDA.
Relationship between LDA, CCA, and RRR
In his answer, @ttnphns made a connection to canonical correlation analysis (CCA). Indeed, LDA can be shown to be equivalent to CCA between $\X$ and $\G$. In addition, CCA between any $\newcommand{\Y}{\mathbf Y}\Y$ and $\X$ can be written as RRR predicting whitened $\Y$ from $\X$. The rest follows from this.
Bibliography
It is hard to say who deserves the credit for what is presented above.
There is a recent conference paper by Cai et al. (2013) On The Equivalent of Low-Rank Regressions and Linear Discriminant Analysis Based Regressions that presents exactly the same proof as above but creates the impression that they invented this approach. This is definitely not the case. Torre wrote a detailed treatment of how most of the common linear multivariate methods can be seen as reduced rank regression, see A Least-Squares Framework for Component Analysis, 2009, and a later book chapter A unification of component analysis methods, 2013; he presents the same argument but does not give any references either. This material is also covered in the textbook Modern Multivariate Statistical Techniques (2008) by Izenman, who introduced RRR back in 1975.
The relationship between LDA and CCA apparently goes back to Bartlett, 1938, Further aspects of the theory of multiple regression -- that's the reference I often encounter (but did not verify). The relationship between CCA and RRR is described in the Izenman, 1975, Reduced-rank regression for the multivariate linear model. So all of these ideas have been around for a while.
