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

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!

Answer 1

[David Goodmanson]

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

[Bruno Luong]

Great David!!

I was on the same track but couldn't figure out the formula.

U = na*u*na' + (1/n)*onzc*onzc';

I was playing with left or right product but did not think about doing in bothside.

[David Goodmanson]

Hi Bruno, yeah, I was feeling the pressure (in a good way) because I figured you were probably on the job.

[Paul]

The OP asked for a "complex unitary" matrix, which suggests the desired answer would have non-zero imaaginary part.

Also, would be interesting to know from OP exactly what is meant by a random matrix. Is there some expectation that the entries of the matrix have some specified joint distribution?

[Bruno Luong]

randc In David code generates complex, presumably is (TBC by David)

randc(n-1) := rand(n-1) + 1i*rand(n-1)

When one deal random generator of elements of a (Lie) group, if the distribution is not specfied it is implicitly understood as uniform wrt the Lie-group Riemannian geometry metric.

The distribution of David code is indeed not uniform if RANDC is as above.

However it becomes uniform with

h = randn(n-1) + 1i*randn(n-1)

[Paul]

I was curious about randc after I saw David's code. The help doesn't say too much. It appears to be an obsolete function. And it returns a real result, at least in 2019a

>> help randc
  randc returns random matrix with entries in [-0.5,0.5]

>> which randc -all
C:\Program Files\MATLAB\R2019a\toolbox\robust\rctobsolete\lmi\randc.m
>> randc(2)

ans =

   3.8794e-01   3.0208e-01
   2.5319e-01  -5.8594e-02

Are you sure you meant randn? I tried using your definition of h and neither the real nor imaginary parts of U(1,1) look uniformly distributed.

[Bruno Luong]

"Are you sure you meant randn"

Yes

"I tried using your definition of h and neither the real nor imaginary parts of U(1,1) look uniformly distributed. "

I never say U(1,1) is uniform.

If you want to generate a vector in r in R^2 such that |r| = 1 (a circle then), do you expect r(1) is uniform on (-1,1)?

[Paul]

Well, I suppose one could want to generate a complex number z with unity norm and angle distributed uniformly on 0-2pi. It doesn't look like that's what's happening here, oor at least I'm not seeing it.

So to be clear: After I apply David's code with h = rand(n-1) + 1i*rand(n-1), what exactly about U will be distributed uniformly?

[David Goodmanson]

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).

[Paul]

David,

Thanks for the update. For sure your method meets the stated requirements of the OP.

[Bruno Luong]

@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.

[Paul]

Thanks for the detailed response.

[David Goodmanson]

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]

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

[Bruno Luong]

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

With this setup

m = randn(n,n) + i*randn(n,n);
% four options
[u ~] = qr(m);

It returns plateau of 2 on my MATLAB (R2020b)

[David Goodmanson]

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]

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]

nor do I.

[Bruno Luong]

Hi @David Goodmanson if you are stll reading. Pusing a little bit further the test, I think RAND with EIG does actually not generate uniform Unitary matrices, as you can see in this graph

obtained from this code (warning it takes relatively long to run)

close all
clearvars
% Generate uniform randomly p Orhogonal/Hermitian matrix in R^n
p = 1e7;
n = 3;

rfun = {@(varargin) randn(varargin{:}), ...
        @(varargin) 2*rand(varargin{:})-1 };

for i = 1:length(rfun)
    randmethod = rfun{i};

    U=randmethod(n,n,p);

    % David EIG
    U = U + permute(conj(U),[2 1 3]);
    for k=1:p
        [U(:,:,k),~]=eig(U(:,:,k));
    end
    % random phase
    rsign = 2*(rand(1,1,p)>0.5)-1;
    rsign = repmat(rsign,[n,1,1]);
    U = U.*rsign;


    % Check for n=2 or 3
    U1=reshape(U(:,1,:),[n p]);

    x = real(U1(1,:));
    y = real(U1(2,:));
    tt = atan2(y,x);

    ax = subplot(1,2,i);
    histogram(ax, tt,100)
    if i == 1
        title(ax,'randn')
    else
        title(ax,'rand')
    end
    drawnow

end