p1<-4 ; p2<-3 ; n<-200
# create a matrix Y with separable covariance
Sig1<-rWishart(1,p1,diag(p1))[,,1]
Sig2<-rWishart(1,p2,diag(p2))[,,1]
Y<-array(rnorm(n*p1*p2),dim=c(n,p1,p2))
Y<-aperm( apply(Y,c(1,3),function(y){ msqrt(Sig1)%*%y } ),c(2,1,3))
Y<-aperm( apply(Y,c(1,2),function(y){ msqrt(Sig2)%*%y } ),c(2,3,1))
# covariance
S<-mcov(Y)
image(S)
plot(S,kronecker(Sig2,Sig1)) ; abline(0,1)
Run the code above in your browser using DataLab