Matrix expression for $\operatorname{vec}(X^{\top}X)$? Let $X$ be an $m$ by $n$ matrix. I would like to find the $n^2$ by $(mn)^2$ matrix $B$ such that
$$
\operatorname{vec}\left(X^{\top}X\right) = B \operatorname{vec}\left(\operatorname{vec}\left(X\right)\left(\operatorname{vec}\left(X\right)\right)^{\top}\right).$$
Here $\operatorname{vec}(\cdot)$ is the vector function that unrolls a matrix into a vector columnwise.
I would guess $B$ involves the Kronecker products, and some $I_m$ and $I_n$ matrices and some 1 vectors, but I am having a hard time with all the indices.
 A: I believe $B$ takes the form
$$
B = \left(I_n \otimes \left(\operatorname{vec}\left(I_m\right)\right)^{\top} \otimes I_n\right)\left(
I_{mn} \otimes K \right),
$$
where $K$ is the 'commutation matrix.' 
That is, $K$ is the matrix such that
$K \operatorname{vec}(X) = \operatorname{vec}\left(X^{\top}\right).$
Letting 
$$
V = \left(I_{mn} \otimes K \right) 
\operatorname{vec}\left(\operatorname{vec}\left(X\right)\left(\operatorname{vec}\left(X\right)\right)^{\top}\right)
= \operatorname{vec}\left(\operatorname{vec}\left(X^{\top}\right)\left(\operatorname{vec}\left(X\right)\right)^{\top}\right),
$$
then the $i,j$ element of $\operatorname{vec}\left(X^{\top}X\right)$ should be:
$$
\operatorname{vec}\left(X^{\top}X\right)_{i,j} = \sum_k X_{k,i} X_{k,j} =
V_{i+mk+mn(k+mj)},
$$
where we index from zero, as suggested by @Chrystomath. The term
$\left(I_n \otimes \left(\operatorname{vec}\left(I_m\right)\right)^{\top} \otimes I_n\right)$
then captures that indexing and the repeated $k$.
I have confirmed this relationship with some R code:
library(matrixcalc)
m <- 5
n <- 3
set.seed(1234)
X <- matrix(rnorm(m*n),nrow=m)
LHS <- vec(t(X) %*% X)
VM  <- vec(vec(X) %*% t(vec(X)))
B   <- (diag(n) %x% t(vec(diag(m))) %x% diag(n)) %*% (diag(m*n) %x% commutation.matrix(m,n))
RHS <- B %*% VM
max(abs(RHS - LHS))

[1] 0

A: Mind boggling indices!
It's better to use the computing convention of running indices $i=0,\ldots,n-1$, $j=0,\ldots,n-1$, $k=0,\ldots,m-1$. 
Let $Y=\mathrm{vec}(\mathrm{vec}(X)\mathrm{vec}(X^T))$. Then $Y_{mn(jm+k)+(im+l)}=x_{il}x_{jk}$.
Similarly, $Z=\mathrm{vec}(X^TX)$ and $Z_{ni+j}=\mathbf{x}_i\cdot\mathbf{x}_j=\sum_{k=0}^{m-1}x_{ik}x_{jk}$.
So, to transfer from $Y$ to $Z$ requires the matrix $$\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}\sum_{k=0}^{m-1}\delta(ni+j,mn(jm+k)+(im+k))$$
That is, the matrix $B$ is that $n^2\times(mn)^2$ matrix with an entry of $1$ in row $ni+j$ and column $mn(jm+k)+(im+k)$ (for the stated ranges of $i,j,k$), and $0$s elsewhere.
I've tested this out on Mathematica:
m=5;n=8;X = Transpose@RandomInteger[{-5, 5}, {m, n}]
Y = Flatten[Transpose[Outer[Times, Flatten[X], Flatten[X]]]]
Z = Flatten@Table[X[[i]].X[[j]], {i, 1, n}, {j, 1, n}]
B = SparseArray@
   Flatten@Table[{n i + j + 1,m n (m j + k) + n i + k + 1} -> 1, {i, 0,
       n-1}, {j, 0, n-1}, {k, 0, m-1}]
Then B.Y matches with Z. (Note that Mathematica does not use the computing convention hence the $+1$ in some of the indices.)
A: So far you've confirmed that your formula works by testing it with various random matrices. That's great, but here is a more theoretical derivation.
Define the quantities
$$\eqalign{
X &\in {\mathbb R}^{m\times n}, \quad &I \in {\mathbb R}^{n\times n} \cr
c_k &= X\varepsilon_k, \quad &r_j = X^Te_j, \quad &x = {\rm vec}(X) \cr
e_j,c_k &\in {\mathbb R}^{m\times 1}, \quad
 &r_j,\varepsilon_k \in {\mathbb R}^{n\times 1}, \quad
 &x \in {\mathbb R}^{mn\times 1} \cr
}$$ where $\{e_j,\varepsilon_k\}$ are the cartesian basis vectors for $\{{\mathbb R}^m,{\mathbb R}^n\}$, $c_k$ is the $k^{th}$ column of $X$, and $r_j^T$ is the $j^{th}$ row. You can also think of $\varepsilon_k$ as the $k^{th}$ column of $I$.
Note that the vectors above can be used to expand various matrix expressions
$$\eqalign{
X &= \sum_{k=1}^n c_k\varepsilon_k^T = \sum_{j=1}^m e_jr_j^T \cr
X^TX &= \sum_{j=1}^m r_jr_j^T \cr
XX^T &= \sum_{k=1}^n c_kc_k^T \cr
}$$
Often I'll just use repeated indices instead of explicit summations, e.g.
$$\eqalign{
{\rm vec}(X^TX) &= {\rm vec}(r_jr_j^T) = r_j\otimes r_j \cr
{\rm vec}(xx^T) &= x\otimes x \cr
}$$
You wish to find a matrix $B$ such that
$$\eqalign{
r_j\otimes r_j &= B\,(x\otimes x)\cr
}$$
If we can find matrices $\{H_j\}$ such that $r_j = H_jx$ then 
$$\eqalign{
B = \sum_{j=1}^m H_j\otimes H_j \cr
}$$ this should be feasible since $x$ contains all the elements of $X$, we just need $H_j$ to select the elements of the $j^{th}$ row. In fact,
$$\eqalign{
H_j &= I\otimes e_j^T \cr
H_jx
 &= (I\otimes e_j^T)\,{\rm vec}(X) \cr
 &= {\rm vec}(e_j^TXI) = {\rm vec}(r_j^T) = r_j \cr
}$$
Therefore the desired coefficient matrix is
$$\eqalign{
B = \sum_{j=1}^m I\otimes e_j^T\otimes I\otimes e_j^T \cr
}$$
Since $I_m = \sum_je_je_j^T\,$ we have
$$\eqalign{
{\rm vec}(I_m)^T &= \sum_{j=1}^m(e_j\otimes e_j)^T
  = \sum_{j=1}^me_j^T\otimes e_j^T \cr
}$$
and your formula for the coefficient matrix can be written as
$$\eqalign{
B = \bigg(\sum_{j=1}^m I\otimes e_j^T\otimes e_j^T\otimes I\bigg)\,(I_{mn}\otimes K_{mn}) \cr
}$$
Your formula requires the commutation matrix to reverse the order of factors $(e_j^T\otimes I)$.
A: Given the matrix $X\in {\mathbb R}^{m\times n}$
Define some vectors in terms of $X$.
$$\eqalign{
 x &= {\rm vec}(X), \quad s = {\rm vec}(X^TX), \quad y = x\otimes x \cr
}$$
Then solve the following linear equation for the $B$ matrix.
$$\eqalign{
 s &= B\,{\rm vec}\Big({\rm vec}(X){\rm vec}(X)^T\Big) = By \cr
 B &= sy^+ + A(I-yy^+) \cr
}$$
where $A$ is an arbitrary matrix the same size as $B$.
Since there is an explicit formula for the pseudoinverse of a vector, i.e.
$$\eqalign{
 y^+ = \frac{y^T}{y^Ty}
}$$
We can also write the result as
$$\eqalign{
 B &= \frac{sy^T + (y^Ty)A-Ayy^T}{y^Ty} \cr
}$$
