# 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.

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

• It should be ${\rm vec}(I_m)$, right? Otherwise dimensions wont work. – Florian Apr 11 at 10:11
• @Florian : oops, you are right. I got it right in the code, but wrong in the math. – shabbychef Apr 12 at 5:21

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.)

• This is helpful, especially the part about zero-indexing. The reason I would like $B$ expressed in matrix form is that my $X$ is a random variable, and I know the covariance of $\operatorname{vec}\left(\operatorname{vec}(X)\left(\operatorname{vec}(X)\right)^{\top}\right)$ and would like to use it to find the covariance of $\operatorname{vec}\left(X^{\top}X\right)$. – shabbychef Apr 7 at 4:27
• The matrix $B$ is given by the formula. I've edited the answer to make it more explicit. Columns and rows start with 0,1,.... – Chrystomath Apr 7 at 8:23

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)$$.

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 }

• I think you wanted $y=\operatorname{vec}\left(X^{\top}X\right)$. In that case, I have confirmed your solution is correct. However, I was looking for a $B$ that works for all $X$. – shabbychef Apr 10 at 5:13
• No, that's the definition of $s$. The $y$ variable is the Kronecker product of ${\rm vec}(X)$ with itself. All of this is in the very first line of my answer. Please note that the solution is valid for all $X$ but also note that it is not unique since it includes the parameter $A$ which is an arbitrary matrix. – greg Apr 10 at 23:46
• ah, you were right. And my comment about working for "all $X$" was imprecise. What I wanted was a $B$ that was independent of $X$, which I will not actually observe. – shabbychef Apr 12 at 5:23