distribution implement for RandomVariate

Question

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.

Answer 1

No idea what your f[x,y] is, so I've just used f[x_,y_] := x.y. This below will work although generating 10^8 samples is going to take a lot longer:

n = 4;
f[x_,y_] := x.y
distribution = MatrixPropertyDistribution[f[x,y],
  Evaluate[{
     y \[Distributed] GaussianUnitaryMatrixDistribution[n], 
     x \[Distributed] MultinormalDistribution[ConstantArray[0, n], DiagonalMatrix[ConstantArray[1, n]]]
  }]
];

sample = RandomVariate[distribution, 10]