# distribution implement for RandomVariate

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.