Least Squares Solution for a Matrix System with Diagonal Matrix Constraint 
Let's say I have the following equation
      $$AX=B$$
  where $A$ is a $8\times 3$ matrix (known), $X$ is a $3\times3$ "diagonal" matrix which represents the coefficients (unknown) and $B$ is a $8\times 3$ matrix (known).

Whenever I solve for $X$ using least squares $X= (A^TA)^{-1} A^T B$, I get a square matrix which is not diagonal.
Any idea how can I force least square solution matrix to be diagonal? 

Decouple your equations to force diagonalization. A simpler problem demonstrates the basics.
Instead of solving
$$
\left(
\begin{array}{cc}
 a_{11} & a_{12} \\
 a_{21} & a_{22} \\
 a_{31} & a_{32}
\end{array}
\right)
\left(
\begin{array}{cc}
 x_{11} & x_{12} \\
 x_{21} & x_{22} 
\end{array}
\right)
%
=
\left(
\begin{array}{cc}
 b_{11} & b_{12} \\
 b_{21} & b_{22} \\
 b_{31} & b_{32}
\end{array}
\right)
$$
solve instead
$$
\left(
\begin{array}{cc}
 a_{11} & a_{12} \\
 a_{21} & a_{22} \\
 a_{31} & a_{32}
\end{array}
\right)
\left(
\begin{array}{cc}
 x_{1} & 0 \\
 0 & x_{2} \\
\end{array}
\right)
%
=
\left(
\begin{array}{cc}
 b_{11} & b_{12} \\
 b_{21} & b_{22} \\
 b_{31} & b_{32}
\end{array}
\right).
$$
This is two separate problems for $k=1,2$:
$$
\left(
  \begin{array}{c}
    a_{1k} \\
    a_{2k} \\
    a_{3k}
  \end{array}
\right)
 x_{k} 
=
\left(
  \begin{array}{c}
    b_{1k} \\
    b_{2k} \\
    b_{3k}
  \end{array}
\right)
$$
The pseudoinverse of a vector $v$ is $\frac{v^{\mathrm{*}}}{\lVert v \rVert}$. For example
$$
\left(
  \begin{array}{c}
    a_{1k} \\
    a_{2k} \\
    a_{3k}
  \end{array}
\right)^{\dagger}
=
\frac{
\left( 
\begin{array}{ccc}
  a_{1k} & a_{2k} & a_{3k}
\end{array}
\right)
}{\lVert a_{k} \rVert}.
$$
The least squares point solution is
$$
  x_{k} = a_{k}^{\dagger} b_{k} = \frac{ \sum_{j=1}^{3}a_{jk}b_{jk} } {\sqrt{ \sum_{j=1}^{3}a_{jk}a_{jk} }}.
$$
