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:

m <- 5
n <- 3
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
  • $\begingroup$ It should be ${\rm vec}(I_m)$, right? Otherwise dimensions wont work. $\endgroup$ – Florian Apr 11 at 10:11
  • $\begingroup$ @Florian : oops, you are right. I got it right in the code, but wrong in the math. $\endgroup$ – 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.)

  • $\begingroup$ 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)$. $\endgroup$ – shabbychef Apr 7 at 4:27
  • $\begingroup$ 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,.... $\endgroup$ – 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 }$$

  • $\begingroup$ 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$. $\endgroup$ – shabbychef Apr 10 at 5:13
  • $\begingroup$ 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. $\endgroup$ – greg Apr 10 at 23:46
  • $\begingroup$ 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. $\endgroup$ – shabbychef Apr 12 at 5:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.