3

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.

Rodrigo de Azevedo
  • 1
tefkah
  • 33
  • 1
    No, because a random $n\times n$ matrix with uniformly distributed entries whose absolute value is at most one will typically have much larger Hilbert-Schmidt norm (that is the square of the trace of $M^*M$) than $1$. – uniquesolution Jul 08 '18 at 20:52
  • 1
    Is it still slow if you vectorize? Generating 10 million 3x3 matrices and summing their squares should be able to be done much faster – Nick Alger Jul 08 '18 at 21:34

1 Answers1

2

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...
Nick Alger
  • 18,844
  • 1
    +1 It should be noticed, however that the resulting distribution would actually not be uniform in $-1,1$ (marginally, for each real and imaginary part). But I guess that this is in the spirit of what is desired. – leonbloy Jul 09 '18 at 00:58
  • 1
    @leonbloy Could you elaborate? I'm not sure I understand. Wouldn't this generate samples from the same distribution as the rejection sampling scheme in the original question? (I think it would, but maybe I'm wrong?) – Nick Alger Jul 09 '18 at 02:38
  • 1
    Yes, this would coincide with the rejection sampling scheme. And neither of them would result in uniformly distributed values for each component. – leonbloy Jul 09 '18 at 02:58
  • 1
    @leonbloy Hmm this is a good point. We will need the question asker to clarify what is wanted – Nick Alger Jul 09 '18 at 03:27
  • Thank you very much! I also saw that you answered my question in a different way first, what made you change it? I don't know if this is exactly as required, my probability theory is not exactly up to snuff, but isn't there a difference between picking uniformly distributed and checking if they're in a sphere different than uniformly picking entries on the sphere? – tefkah Jul 09 '18 at 19:28
  • 1
    @Wibbo There were a few stages. First I just did your method (brute force), except with vectorized code (instead of sequential code). This actually worked well enough for real matrices, but it was still too slow on complex matrices. So I deleted that answer. Next, I did this answer, except without code. I said that one would still need to do rejection sampling, checking that samples from the sphere are in the box. Then when I coded it up I realized that the sphere is actually contained in the box, so the rejection sampling part is unnecessary. That leads to the answer how it is now. – Nick Alger Jul 09 '18 at 19:39
  • 1
    @Wibbo (second part) To make things clear, when I am talking about the uniform distribution on the sphere, I mean the probability distribution that has constant density within the sphere and zero density outside the sphere. Same thing with uniform on the box (constant density in the box, zero outside). Hopefully that is what you want too. Anyways yes picking from the uniform distribution on the box, then rejecting samples that are not in the sphere will yield samples that come from the uniform distribution on the sphere. This follows from the chain rule of probabilty, basically. – Nick Alger Jul 09 '18 at 19:44