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

Question

Michael Lachner 2020-11-21

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/655993-how-to-generate-a-random-complex-unitary-matrix-whose-columns-each-sum-up-to-1

评论： David Goodmanson 2020-12-5

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!

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

David Goodmanson 2020-11-22

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/655993-how-to-generate-a-random-complex-unitary-matrix-whose-columns-each-sum-up-to-1#answer_552183

编辑：David Goodmanson 2020-12-4

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

17 个评论
显示 15更早的评论隐藏 15更早的评论

David Goodmanson 2020-11-23

编辑：David Goodmanson 2020-11-23

Hi Paul / Bruno,

Thanks for pointing this out. In the answer I forgot to replace out randc, which is a homemade function that procuces uniformly distributed numbers on {-1, 1} x {-i, i}, similar to what Bruno had deduced. I modified the answer appropriately.

There has been plenty of discussion in the past about matrices of random numbers. I used the first thing I thought of, a unitary matrix U produced by the eigenvectors of a hermitian matrix M = h+h', with h consisting of uniformly distributed complex random numbers. I'm not making any claim that the resulting U has uniformly distributed elements (which in any case was not part of the question), it just looks pretty random.

The Q from a QR decompostion of M, or U,V from an SVD decompositon of M would be other candidates for U. And M could be uniformly distributed or gaussian distributed, etc. Lots of choices for what might produce the 'most random' unitary matrix.

In the real domain, note that if one uses rand to fill a large matrix, one of its eigenvalues is much, much larger than the rest. That's due to rand's expected value of 1/2. Random numbers on the domain shown above have expected value 0 and avoid that issue. (I know U employs eigenvectors and not eigenvalues but it seemed best to stay away from the whole situation).

Bruno Luong 2020-11-23

编辑：Bruno Luong 2020-11-23

@Paul "It doesn't look like that's what's happening here, oor at least I'm not seeing it."

The reason if if you have a random vector

Z = randn(m) + 1i*randn(m)

The covariance of such random variable is eye(m). Any hermissian transformantion of Z.

T = Q*Z

with Q'*Q = eye(m), (Q in U(n)) also have covariance of eye(m). As normal random (vector) variables is fully characterized by the mean and covariance matrix, T is actually the same random raviable than Z. In other word Z is invariant by unitary transformmation.

That shows that Z is "uniform" on unitary transformation (up to radial scaling).

This is method based on QR (cheaper than EIG) and RANDN I use to generate hermitian matrix. The uniform is checked for n==2 with histogram. You can change RANDN to anything you like and observe that RAND will not generate uniform random matrix,just by checking one vector of it.

I also put EIG method for the record. IMPORTANT it needs an extra random sign swap (phase change in complex) for at least one vector of U with the current algorithm if use with sign symmetric RAND or random

% Generate uniform randomly p Orhogonal/Hermitian matrix in R^n 
p = 1e6;
n = 3;
ComplexFlag = false;
UseQRFlag = true;
if ComplexFlag
    rfun = @(varargin) randn(varargin{:}) + 1i*randn(varargin{:}); % complex
else
    rfun = @(varargin) randn(varargin{:});
    %rfun = @(varargin) 2*rand(varargin{:})-1;  % only for EIG
    %rfun = @(varargin) rand(varargin{:});      % NOT this for any method
end
U=rfun(n,n,p);
if UseQRFlag
    for k=1:p
        [U(:,:,k),~]=qr(U(:,:,k));
    end
    % qr gives fix sign determinant (-1)^p, we need to flip randomly a vector
    rsign = 2*(rand(1,1,p)>0.5)-1;
    rsign = repmat(rsign,[n,1,1]);
    U = U.*rsign;
else
    % David EIG
    U = U + permute(conj(U),[2 1 3]);
    for k=1:p
        [U(:,:,k),~]=eig(U(:,:,k));
    end
    % random phase
    if ComplexFlag
        rsign = exp((2i*pi)*rand(1,1,p));
    else
        rsign = 2*(rand(1,1,p)>0.5)-1;
    end
    rsign = repmat(rsign,[n,1,1]);
    U = U.*rsign;
end
% Check for n=2 or 3
U1=reshape(U(:,1,:),[n p]);
if n==2
    x = real(U1(1,:));
    y = real(U1(2,:));
    tt = atan2(y,x);
    
    close all
    subplot(1,2,1)
    histogram(tt,100)
    subplot(1,2,2)
    plot(x,y,'.')
    axis equal
elseif n==3
    x = real(U1(1,:));
    y = real(U1(2,:));
    z = real(U1(3,:));
    
    [az,el,~] = cart2sph(x,y,z);
    azi = linspace(-pi,pi,65);
    azi(1) = -Inf; azi(end)=Inf;
    eli = asin(linspace(-1,1,33));
    eli(1) = -Inf; eli(end)=Inf;
    [~,i] = histc(az, azi);
    [~,j] = histc(el, eli);
    C = accumarray([i(:) j(:)], 1, [length(azi) length(eli)]-1);
    close all
    subplot(1,2,1)
    bar3(C,'c');
    subplot(1,2,2)
    plot3(x,y,z,'.')
    axis equal  
end

As I say above use RANDN with QR/EIG will generate uniform distribution of U(n) group. This does not imply that U(i,j) is a uniform distribution for any element (i,j) considered alone.

David Goodmanson 2020-12-4

编辑：David Goodmanson 2020-12-4

Hello all,

interesting plots from Bruno, showing uniformity on the sphere for 3x3 real unitary matrices. I treated another aspect here, all the matrix elements of complex U in the complex plane. So consider an nxn complex unitary U, calculated with qr or eig. For each of its elements u, let r = abs(u). Since its columns are normalized to 1, for each column rms(r) = 1/sqrt(n). The same is true for the rms value of all n^2 elements of of U which must be exactly 1/sqrt(n) with no statistical fluctuations.

Plots show that the elements of U form a round blob in the complex plane. For n = 100, the result looks kind of like a globular star cluster:

For larger n such as 2000, you get a solid disc of approximate uniform density with an 'evaporative layer' on the surface. The characteristic radius of the disc scales like 1/sqrt(n), so the density scales like n^3.

It's useful to consider a disc of radius R with uniform density inside, 0 outside, and with rms(r) = 1/sqrt(n). Its radius is R = sqrt(2)/sqrt(n), and the density is rhoR = n^2/(pi*R^2) = n^3/(2*pi). Plot 2 calculates density in relation to rhoR, and is best done with larger values of n such as 2000.

Since an actual distribution of points shows many of them are outside radius R, the density inside R must be larger than rhoR so that the rms value of the whole works stays constant. I didn't use histograms for density, but sorted the r values and did a moving average density plot of rbar vs. N(rbar)/(pi*rbar^2). Here N(rbar) is the number of elements with r <= rbar. If the density is uniform, after some statistical fluctuations near the origin the density should settle down to a horizontal line of height rhoExpt, which it does. Interestingly, the value depends on which method is used, and

rhoExpt/rhoR = approx.  5/3      for qr
rhoExpt/rhoR = approx.  2        for eig

The value does not depend on whether the initial m uses uniform complex or gaussian random complex variables,

n = 2000;
R = sqrt(2)/sqrt(n);
rhoR = n^3/(2*pi);       % rhoR = n^2/(pi*R^2)
% -----------------
% two options
m = (2*rand(n,n)-1) + i*(2*rand(n,n)-1);
%m = randn(n,n) + i*randn(n,n);
% two options
[u ~] = qr(m);
%[u ~] = eig(m+m');
% ------------------
u = u(:);
figure(1)
scatter(real(u),imag(u),1,rand(n^2,1))
map = [repmat([0 0 1],5,1);repmat([1 1 1],20,1);[1 1 0];];
colormap(map)
set(gca,'color','k')
axis('equal')
rbar = sort(abs(u));
movingrho = (1/rhoR)*(1:n^2)'./(pi*rbar.^2); % normalized to rhoR
figure(2)                              
loglog(rbar,movingrho)
grid on
rmsU = rms(u)
1/sqrt(n)    % must equal rmsU

Bruno Luong 2020-12-4

编辑：Bruno Luong 2020-12-4

Hi David,

There is a bug in my code using QR since the phase of the column vectors are not idependent and it shows in the fact that the determinant of u is (-1)^(n-1) (discussed here with Christine Tobler in this thread).

I don't quite follow the argument of density "plateau", but if one introduce a randome phase after qr, it will have a plateau of rhoExpt/rhoR = approx. 2.

I notice that MATLAB EIG implementation introduces a phase shift such that for each eigen vector the component with largest modulus umax becomes real >= 0.

r = rand(10)+1i*rand(10);
[u,~] = eig(r);
[~,loc] = max(abs(u),[],1) % largest modulus
umax = u(sub2ind(size(u),loc,1:size(u,2))) % check if umax is real >= 0

So even it it's less predictable than QR phase of Q, I would also not comfortable of using U as it is after EIG without shuffle the phase if uniformity is desired.

Here is your code modifield with phase shuffle. The RandUGroup function used bellow is from the fex I just upload recently.

n = 2000;
R = sqrt(2)/sqrt(n);
rhoR = n^3/(2*pi);       % rhoR = n^2/(pi*R^2)
% -----------------
% two options
m = (2*rand(n,n)-1) + 1i*(2*rand(n,n)-1);
%m = randn(n,n) + i*randn(n,n);
% four options
%[u ~] = qr(m);
%[u ~] = eig(m+m');
u=RandUGroup(n,'U','UseQRFlag','qr');   % qr with shuffle phase
%u=RandUGroup(n,'U','UseQRFlag',false);  % eig with shuffle phase
% ------------------
u = u(:);
figure(1)
scatter(real(u),imag(u),1,rand(n^2,1))
map = [repmat([0 0 1],5,1);repmat([1 1 1],20,1);[1 1 0];];
colormap(map)
set(gca,'color','k')
axis('equal')
rbar = sort(abs(u));
movingrho = (1/rhoR)*(1:n^2)'./(pi*rbar.^2); % normalized to rhoR
figure(2)                              
loglog(rbar,movingrho)
grid on
rmsU = rms(u)
1/sqrt(n)    % must equal rmsU

David Goodmanson 2020-12-4

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 2020-12-5

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.