I have a scalar function of two random matrices of dimension $ n$ which are in the Gaussian Unitary Ensemble $ $ f(x,y) \in \mathbb R, \quad x^{\dagger}=x \: \operatorname{and} \: y^{\dagger} = y, \:\operatorname{dim}(x)=\operatorname{dim}(y)=n$ $ I would like to generate a sample of various of values of $ f(x,y)$ for example, $ 10^8$ of them. To do this I use the following Mathematica program

`RandomVariate[MatrixPropertyDistribution[f[x, y], {x \[Distributed] GaussianUnitaryMatrixDistribution[n], y \[Distributed] GaussianUnitaryMatrixDistribution[n]}],10^8] `

However, this program runs slow for large dimension $ n$ due to the complicated form of $ f(x,y)$ . So I work in a basis that the matrix $ x$ is diagonalized: $ $ x\rightarrow\operatorname{diag}(a_1,\dots,a_n), \: y \rightarrow U^\dagger y U,$ $ where $ U$ is the unitary matrix that diagonalizes $ x$ . The scalar function $ f$ is invariant under this transformation.

Now $ \{a_1,\dots,a_n\} $ is in a multinormal distribution wherein the covariance matrix is the identity matrix. At the same time, $ y$ is still in the Gaussian Unitary Ensemble. So I would like to generate a sample of values of $ f$ in terms these new random variable. But I find it is hard to implement the distribution. I tried

`distx = Array[x, n] \[Distributed] MultinormalDistribution[ConstantArray[0, n], DiagonalMatrix[ConstantArray[1, n]]]; disty = y \[Distributed] GaussianUnitaryMatrixDistribution[n]; distribution = MatrixPropertyDistribution[f[x, y], Join[distx, {disty}]]; sample = RandomVariate[distribution,10^8]; `

However, I got the following message from Mathematica

Is there any way to solve this problem so that I can obtain the values from the random variables? Thanks very much.