Density Estimation of a Matrix-valued Random Variable? It seems like the density estimation of a multivariate vector-valued random variable has been well studied, but what if one would like to estimate the probability density of a matrix-valued random variable? Would it make sense to perfom a "matrix" form of kernel density estimator, and use a Gaussian kernel function like https://en.wikipedia.org/wiki/Matrix_normal_distribution?
It not, what would be the proper way to do so?

This could make sense.
Matrices are themselves vectors. For $M$ of shape $m\times n$ you can define a vectorization $vec(M) = v_M$ that creates a column vector of dimension $mn$ that consists of all the $n$ columns of $M$ put on top of each other.
So you could just consider matrices as vectors and create Gaussian KDEs as usual.
Now, you referred to this Wikipedia article about "Matrix normal distributions". In there, they prove, that for a given dimension $n\times p$, the matrix random variable $X$ is $\mathcal{MN}_{n\times p}(\mathbf M, \mathbf U, \mathbf V)$ if and only if its vectorization $vec(X)$ is $\mathcal{N}_{np}(vec(\mathbf M), \mathbf V \otimes \mathbf U)$. Now, since the set of covariance matrices, that can be written as $\mathbf V\otimes \mathbf U, V\in\mathbb R^{p\times p}, U\in\mathbb R^{n\times n}$, is a subset of all $np$-dimensional covariance matrices, using distributions from $\mathcal{MN}_{n\times p}$ is effectively imposing a constraint.
If you do want to constrain the possible covariances that way, then your suggestion is the way to go.
In the Wikipedia article, they show that random matrices $X\in\mathbf R^{n\times p}$ that consist of $p$ $n$-dimensional columns that are samples of an $n$-dimensional Gaussian, are $\mathcal{MN}_{n\times p}$. Thus, if you were e.g. dealing only with this type of random matrices, then it would make sense to use a KDE which is a mixture of only $\mathcal{MN}_{n\times p}$ components.
