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? 
 A: Here is the thing. If $X=diag(x_1,x_2,x_3)$ then the elements of first row of $AX$, i.e. $a_{11}x_1, a_{12}x_2, a_{13}x_3$, are equal to $b_{11}$, $b_{12}$ and $b_{13}$ respectively. Hence $x_1=b_{11}/a_{11}$ and $x_2=b_{12}/a_{12}$...further, the second row gives then $x_1=b_{21}/a_{21}$,...and so forth until the 8th row. So in order for $X$ to be a diagonal it must be that the first column of $B$ is a multiple of the first column of $A$, the second column of $B$ a multiple of the second column of $A$ and likewise for the third (that is $b_{.1}=x_1a_{.1}$, $b_{.2}=x_2a_{.2}$,...). I guess your known matrices do not satisfy that.
I guess what you want to do is to minimize $(AX-B)'(AX-B)$ with respect to $X$ under the necessary side conditions (being 'off diagonal elements of $X$ are zero') 
Edit: actually, $X$ being diagonal gives you three equations of the form $b_{.i}=a_{.i}x_i$ for each column $i$. The least square solution is
$$x_i=(a_{.i}'a_{.i})^{-1}a_{.i}'b_{.i} $$ for each $i$.
Now, solving this way means you are treating the columns as separate equations. If they are related by correlated error terms, google for SUR (seemingly unrelated) equations, to find efficient estimators...
A: 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} }}.
$$
A: The solution can be written as
$$\eqalign{
X &= {\rm Diag}\Big((I\odot A^TA)^{-1}\,{\rm diag}(A^TB)\Big) \\
}$$
where $(\odot)$ denotes the elementwise/Hadamard product, Diag() creates a diagonal matrix from a vector argument, and diag() returns the main diagonal of its matrix argument as a vector.
This result is derived as follow.

The independent variable in this problem is the vector $x,\,$ from which a diagonal matrix is defined strictly for notational convenience.
$$\eqalign{
X &= {\rm Diag}(x) \quad\implies\quad  x = {\rm diag}(X) \\
}$$
Given the matrices $(A,B),\,$ the problem is to minimize the function 
$$\eqalign{
\phi &= \tfrac{1}{2}\,\|AX-B\|_F^2 \\
  &= \tfrac{1}{2}\,(AX-B):(AX-B) \\
}$$
where the colon denotes the trace/Frobenius product, i.e. 
$$A:B={\rm Tr}\big(A^TB\big)$$
Calculate the gradient of the function.
$$\eqalign{
d\phi
  &= (AX-B):A\,dX \\
  &= \big(A^TAX-A^TB\big):dX \\
  &= \big(A^TAX-A^TB\big):{\rm Diag}(dx) \\
  &= {\rm diag}\big(A^TAX-A^TB\big):dx \\
\frac{\partial\phi}{\partial x} &= {\rm diag}\big(A^TAX-A^TB\big) \\
}$$
Set the gradient to zero and solve for the optimal vector.
$$\eqalign{
{\rm diag}(A^TAX\,I) &= {\rm diag}(A^TB) \\
(I\odot A^TA)\,x &= {\rm diag}(A^TB) \\
x &= (I\odot A^TA)^{-1}\,{\rm diag}(A^TB) \\
X &= {\rm Diag}\Big((I\odot A^TA)^{-1}\,{\rm diag}(A^TB)\Big) \\
}$$
The preceding makes use of the Hadamard-Diag identity, i.e.
$$\eqalign{
{\rm diag}\Big(A\;{\rm Diag}(x)\;B\Big) &= \big(B^T\odot A\big)\,x \\
}$$
