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

enter image description here

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


1 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],
     y \[Distributed] GaussianUnitaryMatrixDistribution[n], 
     x \[Distributed] MultinormalDistribution[ConstantArray[0, n], DiagonalMatrix[ConstantArray[1, n]]]

sample = RandomVariate[distribution, 10]

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.