for the vec-trick of the Kronecker product, we can write
$$ \left(\mathbf{B}^{\top} \otimes \mathbf{A}\right) \operatorname{vec}(\mathbf{X})=\operatorname{vec}(\mathbf{A} \mathbf{X} \mathbf{B}). $$
I am wondering whether this could be extended to the case of higher-order tensors. Now we are going from a vector to a matrix. Can I reshape the following formula
$$ \left(\mathbf{B}^{\top} \otimes \mathbf{A}\right) \mathbf{X} $$
in something similar to the one above, but where now I have a matrix on the right multiplication? I would like going from a 2D matrix to a 3D matrix in a way. Did you ever encounter this problem? Thanks
The matrix $X$ consists of the columns $\{c_k\},\,$ i.e. $\,X = \big[\matrix{c_1&c_2&\ldots&c_n}\big]$
Consider that each column can be associated with a corresponding matrix $\,c_k={\rm vec}(C_k)$
Reversing the "vec-trick" yields $\;(B^T\otimes A) c_k = {\rm vec}(AC_kB)$
Therefore $$\eqalign{ (B^T\otimes A) X &= (B^T\otimes A)\big[\matrix{c_1&c_2&\ldots&c_n}\big] \\ &= \Big[\matrix{ {\rm vec}(AC_1B) &{\rm vec}(AC_2B) &\ldots &{\rm vec}(AC_nB)}\Big] \\ }$$ which is a partition-wise vectorization of the matrix.
By introducing the standard basis vectors $\{e_k\}$ this can be written more compactly $$\eqalign{ (B^T\otimes A) X &= (B^T\otimes A)\;\sum_{k=1}^nc_ke_k^T \\ &= \big(B^T\otimes A\big)\,c_ke_k^T \qquad&\big({\rm Einstein\,Convention}\big) \\ &= {\rm vec}\big(AC_kB\big)\,e_k^T \\ }$$ One could consider $C_k$ matrix to be the $k^{th}$ component of a third-order tensor $\Gamma$ $$C_k = \Gamma e_k\quad\implies\quad\Gamma = C_ke_k^T$$ The third-order tensor $\Omega$ can be used to define the vec-operation $$\eqalign{ \Omega_{ijk} &= \begin{cases} 1\quad{\rm if}\;i=j+rk-r \\ 0\quad{\rm otherwise} \end{cases} \\ \Omega_{ijk}M_{jk} &= {\rm vec}(M)_i \\ \Omega:M &= {\rm vec}(M) \\ M &\in {\mathbb R}^{r\times c},\quad\; \Omega\in{\mathbb R}^{rc\times r\times c} \\ X &= c_ke_k^T = (\Omega:C_k)e_k^T = \Omega:\Gamma \\ }$$ allowing the solution to be compacted even further $$\eqalign{ (B^T\otimes A) X &= {\rm vec}\big(AC_kB\big)\,e_k^T \\ &= \big(\Omega:AC_kB\big)e_k^T \\ &= \big(\Omega B^T:AC_k\big)e_k^T \\ &= \Omega B^T:A\Gamma \\ }$$ Of course, I don't recommend doing any of this, since these operations are not well-known, or as simple the original Kronecker product.