# 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.

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]