# How to generate a random complex unitary matrix whose columns each sum up to 1

83 views (last 30 days)
Michael Lachner on 21 Nov 2020
Hello everyone,
basically the problem is exactly what is stated in the title. I want to create an unitary (2D) matrix of random complex numbers, such that the elements of each column of this matrix sum up to exactly 1. Is there any way to do this? As long as it's still reasonable, computation speed and/or memory usage are not that important.
Thank you to everyone trying to help!

David Goodmanson on 22 Nov 2020
Edited: David Goodmanson on 4 Dec 2020
MODIFIED to replace previous random function
Hi Michael,
here is one way. It's based on the idea that if the unitary matrix U is nxn, and onz = [1 1 1 1 1 1... ] (length n), then the sum-of-each-column condition is
[1 1 1 1 1 1... ]*U = [1 1 1 1 1 1... ]
so
n = 5;
onz = ones(1,n);
onzc = onz'; % column vector
na = null(onzc');
% construct an (n-1)x(n-1) unitary matrix by employing random numbers
% uniformly distributed on {-1, 1} x {-i, i}
h = 2*(rand(n-1,n-1) + i*rand(n-1,n-1)) -(1+i);
% method 1
[u, ~] = qr(h);
% method 2 (slower)
% [u, ~] = eig(h+h');
% construct the result
U = na*u*na' + (1/n)*onzc*onzc';
% checks, all of these should be small
U'*U -eye(size(U))
U*U' -eye(size(U))
onz*U - onz
onz*U' - onz % U' works too

David Goodmanson on 4 Dec 2020
Hi Bruno,
That's a good point about eig, every eigenvector having one real component, And in your example
r = rand(10)+1i*rand(10);
[u,~] = eig(r);
it's the component with largest absolute value. However for eig(r+r') the real component is usualy not the one with largest absolute value. But just the fact that there is always a real component means that eig needs a phase shuffle.
And I agree that when m is complex gaussian, the plateau value for qr is 2. How did I mess that up?
I probably didn't explain the plateau calculation all that well, The idea is that if you let r = abs(u(:)) and sort the r's, then for every r, all the r's with smaller array index are contained inside a circle of radius r. So for every r, it's a calculation of
rho_average = [number of points < r]/(pi*r^2).
If the density is constant, then rho_average does not change. From the plot there is a large region where this is true (after getting past large fluctuations at small r due to small sample size). Eventually you get to wider scattered points, not uniformly distributed in radius, and rho_average drops down.
I don't agree that [ uniform random variable and qr ] followed by a phase shuffle will raise the plateau value from 1.6 to 2, because phase shuffle has no effect on the r values.
Bruno Luong on 5 Dec 2020
Yes I was wrong on the plateau value of qr since I though the RAND/RANDN do not have the effect, so I run my RandUnitary code and the only thing that add is phase shuffle.
Indeed the drop of plateau value is due to RAND as we both observe now. I knew RAND is not a good input at least for QR.
I'm still stumped that EIG can be somewhat pass few uniformity tests when using with RAND. I admit that I still don't fully understand it.
David Goodmanson on 5 Dec 2020
nor do I.