# Matrix Variate Normal Distributions with MixMatrix

#### Geoffrey Thompson

#### 2019-11-14

    library(MixMatrix)

## Matrix Variate Distributions

Suppose there is a two-dimensional grid of observations: each observation is randomly distributed somehow, perhaps each spot on the grid has its own mean, there’s a covariance relation across the rows that doesn’t vary across the columns and a covariance relation across the columns that doesn’t vary across the rows. One case where this might come up is a multivariate time series: at each point in time we are dealing with a multivariate observation (perhaps normally distributed or not) and then there is a covariance structure across time to account for. We may be interested in drawing from such distributions, calculating densities, or estimating parameters based on observations. This package presents some functions for doing so in the case of matrix variate normal distributions. For more details about matrix variate distributions in general and matrix variate normal distributions in particular, see Gupta and Nagar (1999), whose results I rely on for much of this presentation.

## Matrix Variate Normal Distributions

A random matrix \(\mathbf{X}\) \((p \times n)\) with a matrix variate normal distribution is parameterized by a mean matrix \(\mathbf{M}\) \((p \times n)\) and covariance matrices \(\mathbf{U}\) \((p \times p)\) (determining covariance among rows) and \(\mathbf{V}\) \((n \times n)\) (determining covariance among columns) and has the following probability density function (Wikipedia (2018)): \[p(\mathbf{X}\mid\mathbf{M}, \mathbf{U}, \mathbf{V}) = \frac{\exp\left( -\frac{1}{2} \, \mathrm{tr}\left[ \mathbf{V}^{-1} (\mathbf{X} - \mathbf{M})^{T} \mathbf{U}^{-1} (\mathbf{X} - \mathbf{M}) \right] \right)}{(2\pi)^{np/2} |\mathbf{V}|^{n/2} |\mathbf{U}|^{p/2}} \] Here is a useful fact for drawing observations from a matrix variate normal distribution: suppose \(\mathbf{X}\) is distributed as a matrix variate normal distribution with mean matrix \(\mathbf{0}\) and covariance matrices \(\mathbf{I_n}\) and \(\mathbf{I_p}\). Then for a mean matrix \(\mathbf{M}\) and linear transformations \(\mathbf{L}\) and \(\mathbf{R}\) of appropriate dimensions, if \[Y = \mathbf{M} + \mathbf{L}\mathbf{X}\mathbf{R}\] then \(\mathbf{Y}\) is matrix variate normally distributed with parameters \(\mathbf{M}, \mathbf{LL}^T, \mathbf{R}^T\mathbf{R}\).

Matrix variate random variables can then be generated by specifying \(\mathbf{M}\) and \(\mathbf{L}\) and \(\mathbf{R}\) or by specifying \(\mathbf{U}\) and \(\mathbf{V}\). If the covariance matrices are provided, then Cholesky decomposition is performed to create the \(\mathbf{L}\) and \(\mathbf{R}\) matrices. If \(\mathbf{L}\) and \(\mathbf{R}\) matrices are provided, they are used to create \(\mathbf{U}\) and \(\mathbf{V}\) matrices which are then decomposed. If these are not specified, the variances are presumed to be identity matrices by default and the means are presumed to be zero. If \(\mathbf{L}\) or \(\mathbf{R}\) are not full rank linear transformations, degenerate matrix normal distributions can be produced.

    set.seed(20180203)
    x <- matrix(rnorm(6),nrow=3)
    x
    #>             [,1]      [,2]
    #> [1,]  0.04234611 1.7072872
    #> [2,] -1.05186238 0.2182058
    #> [3,]  0.36199161 0.6759820
    set.seed(20180203)
    y <- rmatrixnorm(n = 1, mean=matrix(rep(0,6),nrow=3))
    y
    #>             [,1]      [,2]
    #> [1,]  0.04234611 1.7072872
    #> [2,] -1.05186238 0.2182058
    #> [3,]  0.36199161 0.6759820
    U <- 5 * diag(3) + 1
    V <- matrix(c(2,0,0,.1),nrow=2)
    mu = matrix(1:6,nrow=3)
    set.seed(20180203)
    z <- rmatrixnorm(n = 1, mean=mu,U=U,V=V)
    z
    #>           [,1]     [,2]
    #> [1,]  1.146691 5.322459
    #> [2,] -1.568345 5.387067
    #> [3,]  3.734947 6.755212
    mu + t(chol(U)) %*% y %*% chol(V)
    #>           [,1]     [,2]
    #> [1,]  1.146691 5.322459
    #> [2,] -1.568345 5.387067
    #> [3,]  3.734947 6.755212

When \(n=1\), by default this returns a matrix. For \(n>1\), by default, it will return a three-dimensional array indexed by the last coordinate (array convention similar to `rWishart`).

    set.seed(20180203)
    x <- rmatrixnorm(n = 1, mean=matrix(rep(0,6),nrow=3))
    x
    #>             [,1]      [,2]
    #> [1,]  0.04234611 1.7072872
    #> [2,] -1.05186238 0.2182058
    #> [3,]  0.36199161 0.6759820
    set.seed(20180203)
    y <- rmatrixnorm(n = 100, mean=matrix(rep(0,6),nrow=3),list = TRUE)
    y[[1]]
    #>             [,1]      [,2]
    #> [1,]  0.04234611 1.7072872
    #> [2,] -1.05186238 0.2182058
    #> [3,]  0.36199161 0.6759820
    set.seed(20180203)
    z <- rmatrixnorm(n = 100, mean=matrix(rep(0,6),nrow=3),array = TRUE)
    z[ , , 1]
    #>             [,1]      [,2]
    #> [1,]  0.04234611 1.7072872
    #> [2,] -1.05186238 0.2182058
    #> [3,]  0.36199161 0.6759820

Densities are computed in much the same way as for other distributions. The defaults for the mean parameters are 0 and identity for the variance matrices. The `log` parameter sets whether to return the density on the log scale - by default this is `FALSE`.

    set.seed(20180202)
    A = rmatrixnorm(n=1,mean=matrix(c(100,0,-100,0,25,-1000),nrow=2),
    L=matrix(c(2,1,0,.1),nrow=2))
    dmatrixnorm(A,mean=matrix(c(100,0,-100,0,25,-1000),nrow=2),
      L=matrix(c(2,1,0,.1),nrow=2),log=TRUE )
    #> [1] -4.36614

## Parameter estimation for matrix variate normal distributions

Maximum likelihood estimation is possible for the matrix variate normal distribution (see Dutilleul (1999) or Glanz and Carvalho (2013) for details). The pair of variance matrices are only identifiable up to a constant, so this function sets the first element of \(\mathbf{U}\) to \(1\). There is a closed-form solution for the mean matrix and the \(\mathbf{U}\) and \(\mathbf{V}\) have an iterative solution.

    set.seed(20180202)
    A = rmatrixnorm(n=100,mean=matrix(c(100,0,-100,0,25,-1000),nrow=2),
       L=matrix(c(2,1,0,.1),nrow=2),list=TRUE)
    results=MLmatrixnorm(A)
    print(results)
    #> $mean
    #>            [,1]         [,2]      [,3]
    #> [1,] 99.7692446 -100.2587675   25.0921
    #> [2,] -0.1010964   -0.1537989 -999.9448
    #>
    #> $U
    #>           [,1]      [,2]
    #> [1,] 1.0000000 0.5011833
    #> [2,] 0.5011833 0.2542832
    #>
    #> $V
    #>             [,1]        [,2]         [,3]
    #> [1,] 1.000000000  0.08886027  0.003307182
    #> [2,] 0.088860275  0.99216701 -0.048960852
    #> [3,] 0.003307182 -0.04896085  0.808693731
    #>
    #> $var
    #> [1] 3.984766
    #>
    #> $iter
    #> [1] 5
    #>
    #> $tol
    #> [1] 9.830401e-11
    #>
    #> $logLik
    #> [1] -415.1583 -376.4587 -376.4574 -376.4574 -376.4574
    #>
    #> $convergence
    #> [1] TRUE
    #>
    #> $call
    #> MLmatrixnorm(data = A)

There are two restrictions possible for the mean matrices: `row.mean = TRUE` will force a common mean within a row and `col.mean = TRUE` will force a common mean within a column. Setting both will ensure a constant mean for the entire system. Restrictions on \(\mathbf{U}\) and \(\mathbf{V}\) are possible with `row.variance` and `col.variance` commands.

    results.fixrows = MLmatrixnorm(A, row.mean = TRUE, max.iter = 5)
    #> Warning in MLmatrixnorm(A, row.mean = TRUE, max.iter = 5): Failed to
    #> converge
    print(results.fixrows)
    #> $mean
    #>            [,1]       [,2]       [,3]
    #> [1,]    8.20086    8.20086    8.20086
    #> [2,] -333.39989 -333.39989 -333.39989
    #>
    #> $U
    #>           [,1]      [,2]
    #> [1,] 1.0000000 0.4732438
    #> [2,] 0.4732438 0.2752951
    #>
    #> $V
    #>           [,1]      [,2]      [,3]
    #> [1,]  1.000000  1.287832 -2.287827
    #> [2,]  1.287832  1.717146 -3.004973
    #> [3,] -2.287827 -3.004973  5.292803
    #>
    #> $var
    #> [1] 441310.5
    #>
    #> $iter
    #> [1] 5
    #>
    #> $tol
    #> [1] 2.815079e+12
    #>
    #> $logLik
    #> [1] -3125.019 -3093.986 -3062.230 -3030.199 -2998.804
    #>
    #> $convergence
    #> [1] FALSE
    #>
    #> $call
    #> MLmatrixnorm(data = A, row.mean = TRUE, max.iter = 5)
    # this failure is expected with misspecification! The number of iterations is also
    # fixed to be low so the vignette compiles quickly.

Currently the options for variance restrictions are to impose an AR(1) structure by providing the `AR(1)` option, a compound symmetry structure by providing the `CS` option, to impose a correlation matrix structure by specifying `correlation` or `corr`, or to impose an identical and independent structure by specifying `Independent` or `I`. This works by using `uniroot` to find the appropriate \(\rho\) which sets the derivative of the log-likelihood to zero for the `AR` and `CS` options - it is not fast but if this is the true structure it will be better in some sense than an unstructured variance matrix. The \(\rho\) parameter should be \(>0\) and is forced to be non-negative. If the data behaves incompatibly with those restrictions, the function will provide a warning and exit with the current model fit.

    # tolerance and maximum iterations are limited here to make the vignette compile faster

    B = rmatrixnorm(n = 50, mean = matrix(1:15,nrow = 5),
                    U = 3*stats::toeplitz(.6^(1:5)))
    MLmatrixnorm(B, row.variance = "AR(1)", tol = 1e-5)
    #> $mean
    #>           [,1]     [,2]     [,3]
    #> [1,] 0.9879945 5.976905 10.88718
    #> [2,] 1.8732469 6.852212 12.18793
    #> [3,] 2.9224880 7.939933 13.30729
    #> [4,] 4.2194033 9.036333 14.06759
    #> [5,] 5.2340302 9.857787 14.63149
    #>
    #> $U
    #>           [,1]      [,2]      [,3]      [,4]      [,5]
    #> [1,] 1.0000000 0.5839879 0.3410418 0.1991643 0.1163095
    #> [2,] 0.5839879 1.0000000 0.5839879 0.3410418 0.1991643
    #> [3,] 0.3410418 0.5839879 1.0000000 0.5839879 0.3410418
    #> [4,] 0.1991643 0.3410418 0.5839879 1.0000000 0.5839879
    #> [5,] 0.1163095 0.1991643 0.3410418 0.5839879 1.0000000
    #>
    #> $V
    #>              [,1]         [,2]        [,3]
    #> [1,] 1.0000000000 0.0001430499 0.005716386
    #> [2,] 0.0001430499 0.8494982020 0.035590626
    #> [3,] 0.0057163858 0.0355906255 0.979331742
    #>
    #> $var
    #> [1] 1.853556
    #>
    #> $iter
    #> [1] 12
    #>
    #> $tol
    #> [1] 8.350665e-06
    #>
    #> $logLik
    #>  [1] -1156.190 -1149.695 -1148.326 -1147.869 -1147.711 -1147.650 -1147.627
    #>  [8] -1147.617 -1147.614 -1147.612 -1147.611 -1147.611
    #>
    #> $convergence
    #> [1] TRUE
    #>
    #> $call
    #> MLmatrixnorm(data = B, row.variance = "AR(1)", tol = 1e-05)
    MLmatrixnorm(B, tol = 1e-5)
    #> $mean
    #>           [,1]     [,2]     [,3]
    #> [1,] 0.9879945 5.976905 10.88718
    #> [2,] 1.8732469 6.852212 12.18793
    #> [3,] 2.9224880 7.939933 13.30729
    #> [4,] 4.2194033 9.036333 14.06759
    #> [5,] 5.2340302 9.857787 14.63149
    #>
    #> $U
    #>           [,1]      [,2]      [,3]      [,4]      [,5]
    #> [1,] 1.0000000 0.5088212 0.3116836 0.2418073 0.1278731
    #> [2,] 0.5088212 0.9878158 0.6383707 0.3910373 0.3114872
    #> [3,] 0.3116836 0.6383707 1.1270834 0.6297656 0.4495454
    #> [4,] 0.2418073 0.3910373 0.6297656 1.0046405 0.6505338
    #> [5,] 0.1278731 0.3114872 0.4495454 0.6505338 1.0516576
    #>
    #> $V
    #>             [,1]        [,2]        [,3]
    #> [1,] 1.000000000 0.005092835 0.009893861
    #> [2,] 0.005092835 0.836800036 0.034013446
    #> [3,] 0.009893861 0.034013446 0.972453799
    #>
    #> $var
    #> [1] 1.8111
    #>
    #> $iter
    #> [1] 4
    #>
    #> $tol
    #> [1] 2.457002e-06
    #>
    #> $logLik
    #> [1] -1145.534 -1143.775 -1143.767 -1143.766
    #>
    #> $convergence
    #> [1] TRUE
    #>
    #> $call
    #> MLmatrixnorm(data = B, tol = 1e-05)

In order to fit parameters for an \(n \times p\) matrix variate normal distribution, at least \(\max (n/p,p/n) + 1\) observations are needed. Maximum likelihood estimates are unique if there are more than \(\max(p,n)\) observations.

The function expects data in a format produced by the `rmatrixnorm` function. The results are only the parameter estimates and do not include measures of uncertainty.

## Session Information

This vignette was compiled with `rmarkdown`.

    sessionInfo()
    #> R version 3.6.1 (2019-07-05)
    #> Platform: x86_64-redhat-linux-gnu (64-bit)
    #> Running under: Fedora 30 (Thirty)
    #>
    #> Matrix products: default
    #> BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so
    #>
    #> locale:
    #>  [1] LC_CTYPE=en_US.utf8       LC_NUMERIC=C
    #>  [3] LC_TIME=en_US.utf8        LC_COLLATE=C
    #>  [5] LC_MONETARY=en_US.utf8    LC_MESSAGES=en_US.utf8
    #>  [7] LC_PAPER=en_US.utf8       LC_NAME=C
    #>  [9] LC_ADDRESS=C              LC_TELEPHONE=C
    #> [11] LC_MEASUREMENT=en_US.utf8 LC_IDENTIFICATION=C
    #>
    #> attached base packages:
    #> [1] stats     graphics  grDevices utils     datasets  methods   base
    #>
    #> other attached packages:
    #> [1] magrittr_1.5    dplyr_0.8.3     ggplot2_3.2.1   MixMatrix_0.2.4
    #>
    #> loaded via a namespace (and not attached):
    #>  [1] Rcpp_1.0.3        knitr_1.25        tidyselect_0.2.5
    #>  [4] munsell_0.5.0     colorspace_1.4-1  R6_2.4.0
    #>  [7] rlang_0.4.0       highr_0.8         stringr_1.4.0
    #> [10] tools_3.6.1       grid_3.6.1        gtable_0.3.0
    #> [13] xfun_0.10         withr_2.1.2       htmltools_0.4.0
    #> [16] assertthat_0.2.1  yaml_2.2.0        lazyeval_0.2.2
    #> [19] digest_0.6.22     CholWishart_1.1.0 tibble_2.1.3
    #> [22] crayon_1.3.4      purrr_0.3.3       glue_1.3.1
    #> [25] evaluate_0.14     rmarkdown_1.16    labeling_0.3
    #> [28] stringi_1.4.3     compiler_3.6.1    pillar_1.4.2
    #> [31] scales_1.0.0      pkgconfig_2.0.3

## All the code for easy copying

    knitr::opts_chunk$set(
      collapse = TRUE,
      comment = "#>"
    )
    set.seed(20180203)
    x <- matrix(rnorm(6),nrow=3)
    x
    set.seed(20180203)
    y <- rmatrixnorm(n = 1, mean=matrix(rep(0,6),nrow=3))
    y
    U <- 5 * diag(3) + 1
    V <- matrix(c(2,0,0,.1),nrow=2)
    mu = matrix(1:6,nrow=3)
    set.seed(20180203)
    z <- rmatrixnorm(n = 1, mean=mu,U=U,V=V)
    z
    mu + t(chol(U)) %*% y %*% chol(V)

    set.seed(20180203)
    x <- rmatrixnorm(n = 1, mean=matrix(rep(0,6),nrow=3))
    x
    set.seed(20180203)
    y <- rmatrixnorm(n = 100, mean=matrix(rep(0,6),nrow=3),list = TRUE)
    y[[1]]
    set.seed(20180203)
    z <- rmatrixnorm(n = 100, mean=matrix(rep(0,6),nrow=3),array = TRUE)
    z[ , , 1]

    set.seed(20180202)
    A = rmatrixnorm(n=1,mean=matrix(c(100,0,-100,0,25,-1000),nrow=2),
    L=matrix(c(2,1,0,.1),nrow=2))
    dmatrixnorm(A,mean=matrix(c(100,0,-100,0,25,-1000),nrow=2),
      L=matrix(c(2,1,0,.1),nrow=2),log=TRUE )

    set.seed(20180202)
    A = rmatrixnorm(n=100,mean=matrix(c(100,0,-100,0,25,-1000),nrow=2),
       L=matrix(c(2,1,0,.1),nrow=2),list=TRUE)
    results=MLmatrixnorm(A)
    print(results)
    results.fixrows = MLmatrixnorm(A, row.mean = TRUE, max.iter = 5)
    print(results.fixrows)
    # this failure is expected with misspecification! The number of iterations is also
    # fixed to be low so the vignette compiles quickly.
    # tolerance and maximum iterations are limited here to make the vignette compile faster

    B = rmatrixnorm(n = 50, mean = matrix(1:15,nrow = 5),
                    U = 3*stats::toeplitz(.6^(1:5)))
    MLmatrixnorm(B, row.variance = "AR(1)", tol = 1e-5)
    MLmatrixnorm(B, tol = 1e-5)

    sessionInfo()

## References

Dutilleul, Pierre. 1999. “The Mle Algorithm for the Matrix Normal Distribution.” Journal of Statistical Computation and Simulation 64 (2). Taylor & Francis: 105–23.

Glanz, H., and L. Carvalho. 2013. “An Expectation-Maximization Algorithm for the Matrix Normal Distribution.” ArXiv E-Prints, September.

Gupta, Arjun K, and Daya K Nagar. 1999. Matrix Variate Distributions. Vol. 104. CRC Press.

Wikipedia. 2018. “Matrix normal distribution — Wikipedia, the Free Encyclopedia.” http://en.wikipedia.org/w/index.php?title=Matrix%20normal%20distribution&oldid=757552257.