Closed pohwanwuzan closed 2 months ago
I'm not a maintainer of this package, but I have recently been studying the same problem in order to sample random unitaries for my own package, and I believe I understand the issue.
First of all, in order to correct the sign of the matrix, you need to transform it from a bunch of Householder reflections to a regular matrix, but this transformation alone is O(d^3), which kind of defeats the purpose of Stewart's algorithm in offering a O(d^2) method. In his paper he proposes storing the unitary as the set of householder reflections together with the diagonal matrix to correct the signs. Which of course works, but some type other than QRPackedQ
would be needed. Nevertheless, Stewart's algorithm is still much faster than doing the full QR decomposition, so it's still worth using it even without going for exotic types.
Secondly, I don't think what you are doing is correct. You can't correct the sign via a scalar. And furthermore Julia's qr
makes the diagonal of R
real even for complex matrices, so there's no need to use complex phases, +-1 suffices.
Note this package doesn’t really have an active maintainer so if you wish to contribute and eventually take over as maintainer it would be really appreciated
I appreciate the vote of confidence, but I'm afraid this package needs more love than I can give. I can fix the implementation of Stewart's algorithm, but I have neither the time nor the knowledge needed for something deeper.
Maybe it's a good idea to advertise the package in JuliaCon? There must be plenty of people familiar with random matrix theory there, and I don't think this is a case where we need to be afraid of Jia Tan.
Hi there,
First of all thanks for this nice package. I had a problem when I tried to use the function
Stewart
to generate unitary matrices in my simulation though. The function rests in the fileHaarMeasure.jl
, so I assume that it's intended to produce unitary matrices distributed with Haar measure. According to [1], the normalized eigenvalue density is $\rho(\theta) = \frac{1}{2 \pi}$, where $\theta$ is the argument of the eigenvalues of the generated matrices.I tried to reproduce Figure 2 in [1] with the following script:
The resulting density does not seem to be constant; instead it looks like the wrong distribution in Figure 1 in [1]:
Looking into the source code, it seems that the function calls
larfg
from LAPACK to generate Householder reflectors, but the reflectors themselves are not rotated as suggested in [1] (specifically equation (7.22)). To test out I just made a simple change to the source code:Then I modified the script and carried out the small experiment again:
The density now appears constant:
Perhaps I'm missing something here, and it would be great to hear your opinion. Cheers!
References
[1] Mezzadri, F., 2006. How to generate random matrices from the classical compact groups. arXiv preprint math-ph/0609050.