Matrix expression for $\operatorname{vec}(X^{\top}X)$?

Question

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.

Answer 1

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

Answer 2

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

Answer 3

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 }$$

Answer 4

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