The \(r \times c\) matrices \({\bf{H}}{}_{i,j}\) constructed
by the function H.matrices are combined using direct product
to generate the commutation product with the formula \({{\bf{K}}_{r,c}} = \sum\limits_{i = 1}^r {\sum\limits_{j = 1}^c {\left( {{{\bf{H}}_{i,j}} \otimes {{{\bf{H'}}}_{i,j}}} \right)} }\)
References
Magnus, J. R. and H. Neudecker (1979). The commutation matrix: some properties and applications,
The Annals of Statistics, 7(2), 381-394.
Magnus, J. R. and H. Neudecker (1999) Matrix Differential Calculus with Applications in Statistics and Econometrics,
Second Edition, John Wiley.