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.