Sample variance matrix normal distribution in R

I'm trying to perform a multinomial logistic regression in R employing the Metropolis-Hastings algorithm, considering a Matrix Normal Distribution as proposal. I'm using the function rmatrixnorm() included in the package LaplacesDemon in order to sample from the proposal distribution. I followed this strategy since I need a vector of parameter $\underline{\beta_{k}}$, with $k=1,\dots,K$ (number of classes involved in the classification).

At the end of the Monte-Carlo iterations, my procedure retrieves the sample mean and the sample covariance of the posterior distribution. My question is: is there any pre-built function that allows to obtain the sample estimations of the matrices $U$ and $V$ of the Matrix Normal distribution?

Here is my code:

  library(LaplacesDemon)
  multinomial_logit - function(X,y,n_classes,mu,tau2,k,R=1e4,seed=12345){
  set.seed(seed)
  
  #Step 1: Same prior for each parameter vector.
  n = nrow(X)
  m = ncol(X)
  La_rows = sqrt(tau2)*diag(n_classes)
  La_cols = sqrt(tau2)*diag(m)
  
  # Setting-up the matrices U, V dividing them by the precision parameter k
  La_rows = La_rows/k
  La_cols = La_cols/k
  acc = 0
  
  be = matrix(0,n_classes,m)
  bes = matrix(0,n_classes,m)
  
  histMatrix = array(dim=c(R,n_classes,m))
  
  for(h in 1:R){
    
    #Sampling from the proposal
    bes = rmatrixnorm(be,La_rows,La_cols)

    #Probabilities estimation (Softmax Function)
    ps = exp(X%*%t(bes))/(sum(exp(X%*%t(bes))))
    p = exp(X%*%t(be))/(sum(exp(X%*%t(be))))
    
    lnum = dmatrixnorm(bes,mu,La_rows,La_cols,log=TRUE)+dcat(y,p=ps,log=TRUE)
    lden = dmatrixnorm(be,mu,La_rows,La_cols,log=TRUE)+dcat(y,p=p,log=TRUE)
    
    #Testing the acceptance condition of the M-H algorithm
    al = min(1,exp(lnum-lden))
    r = runif(1)
    if(r=al){
      acc = acc+1
      be = bes
    }
    
    histMatrix[h,,] = be
  }
  
  print(dim(histMatrix))
  
  # output (WIP)
  acc_rate = acc/R
  mut = apply(histMatrix[-(1:(round(R/4)+1)),,], c(2,3), mean)
  #Here I want to estimate the matrices U and V from the sample I have generated!
  tau2t = apply(histMatrix[-(1:(round(R/4)+1)),,], c(2,3), var) #WRONG!!!
  return(list(approximated_mean=as.vector(mut),posterior_variance=tau2t,acc_rate=acc_rate))
}
```