Is there a way to quickly generate a random $n\times n$ matrix $M$ with uniformly distributed entries with $\mbox{Tr} (M^\dagger M) \leq 1$?

Question

I am trying to reproduce the results of a paper where they generate $10^6$ $3\times 3$ matrices $M$ with complex elements both real and imaginary uniformly distributed from $-1$ to $1$. They impose another condition, namely that $\mbox{Tr} (M^\dagger M) \leq 1$. I have tried to generate matrices that fulfill this condition by first generating one and checking whether it does, if not I generate a new one etc.

This however takes extremely long, with one matrix taking as long as ten million tries to find, which on my not so stellar PC is upwards of a couple minutes. I don't have months to run this program, so I was wondering whether there are some constraints I can put on the randomly generated entries that would make sure that I would get a matrix that satisfies this condition.

After generating a few I did notice that the norm of the entries usually didn't exceed $0.5$, however this was not always the case. I want to bias my sampling as little as possible, so if there is something to be said about the norm I need to know what this could be.

I tried setting the norm of the elements to be less than $0.5$, $0.6$, $0.7$, etc, and while that did speed it up immensely I have yet to find a maximum.

A solution I did think of is that since $\mbox{Tr} (M^\dagger M) \leq 1$, $\sum_{i,j} |M_{ij}|^2 \leq 1$ and thus the average of the norm of every element has to be smaller than $\frac{1}{3}$, but I don't know how to use this without biasing the trace towards 1.

Answer 1

If I understand this question correctly, this reduces to generating samples from the unit ball in $\mathbb{R}^{2k^2}$, where the matrices are $k\times k$. For concreteness we take $k=3$. However, there is some confusion about what is wanted (see comments below), perhaps the question asker could clarify.

First, reduce the problem from matrices to vectors. Since $\text{trace}(M^*M)$ is just the sum of squares of the absolute values of the entries of $M$, picking $3\times 3$ matrices you want reduces to the problem of picking uniformly random points from $\mathbb{C}^9$ that are both within the box $([-1,1] + i[-1,1])^9$, and within the unit ball in $\mathbb{C}^9$. Actually the constraint that the points are within the box is superflous, since any point in the unit ball will also be in the box.

Next reduce the problem from complex numbers in $\mathbb{C}^9$ to real numbers in $\mathbb{R}^{18}$. You must pick uniformly random points from $\mathbb{R}^{18}$ that are within the unit ball in $\mathbb{R}^{18}$.

The rejection sampling scheme you are using picks random points from within the box, then checks if they are within the ball. Instead, you could pick random points from within the ball directly. To do this, you can use the scheme described by Nate Eldridge here: https://math.stackexchange.com/a/87238/3060

Basically, you pick a random point on the surface of the unit sphere by picking a gaussian random vector and normalizing it. Then you pick a random radius $r$ by picking a uniformly random point $u$ in $[0,1]$, and scaling it: $r=u^{1/d}$.

Here is vectorized MATLAB/Octave code to solve the problem using the method outlined above. It generates $10^6$ sample matrices on my computer in less than a second.

ns = 1e6; % number of samples
n=3; % Matrices will be n-by-n

normalize = @(V) V./repmat(sqrt(sum(abs(V).^2, 1)), size(V,1),1);
randball = @(dim, num_samples) repmat(rand(1,num_samples).^(1/dim), dim, 1) ...
    .* normalize(randn(dim, num_samples));

S = reshape(randball(2*n*n, ns), 2,n,n,ns);
matrices = squeeze(S(1,:,:,:) + 1j*S(2,:,:,:));

first_matrix = matrices(:,:,1) % for example...