Incompability of Krochener products between `rmatrixnormal()` and `rmatnorm()`

[phargarten2]

Greetings! I am the author of matrixNormal package, and I also wrote a function to generate random numbers for the matrix Normal distribution (rmatnorm). I used mvtnorm::rmvnorm() to sample, but the critical question is how the covariance matrix, Sigma, is defined. Unfortunately, we both use the notation U and V but with different meanings so I will use the description (row/column) instead. I am writing to ask about the use of Krochener product of column-row instead of row-column. A user recently asked me about this: #https://github.com/phargarten2/matrixNormal/issues/1. I doubled check this cited reference in the help file (rmatnorm); indeed, Kronecker (row, column) appears to be correct.

Iranmanesh, Anis, M. Arashi, and S. M. M. Tabatabaey On Conditional Applications of Matrix Variate Normal Distribution. Iranian Journal of Mathematical Sciences and Informatics 5, no. 2. (November 1, 2010): 33-43. < https://doi.org/10.7508/ijmsi.2010.02.004 >

Obviously, there is an error between our two functions, and I may be in the wrong. Could you give some insight in answering the user's question specifically and in constructing the rmatrixnormal() function generally? Thanks!

[stla]

Hello,

I did this package a couple of years ago and I don't remember a lot of things. I'm not sure I get your question.

I've checked that E[X X'] = tr(V) U when the mean matrix is null (this formula is given on Wikipedia):

library(matrixsampling)

m <- 3
p <- 2
M <- matrix(0, m, p)
U <- toeplitz(m:1)
V <- toeplitz(p:1)
MNsims <- rmatrixnormal(50000, M, U, V)

X <- matrix(0, nrow = m, ncol = m)
for(i in 1:50000){
  X <- X + tcrossprod(MNsims[, , i])
}
X/50000
sum(diag(V)) * U

This gives:

> X/50000
          [,1]      [,2]      [,3]
[1,] 11.970516  7.990219  4.016199
[2,]  7.990219 11.980130  7.983620
[3,]  4.016199  7.983620 11.943641
> sum(diag(V)) * U
     [,1] [,2] [,3]
[1,]   12    8    4
[2,]    8   12    8
[3,]    4    8   12

So I believe my implementation is correct.

[phargarten2]

Thank you. I hope that you were not offended by my question. Yes, I was wondering about the reference of your function. For your formula, did you mean that *E[X'X] = tr(V) U* ? Thanks to @prockenschaub, I realize that the Iranmanesh reference may not be correct. The evidence seems to indicate that it is the column row variance instead, as Pocuca et al. (2019) suggests.