이 질문을 팔로우합니다.
- 팔로우하는 게시물 피드에서 업데이트를 확인할 수 있습니다.
- 정보 수신 기본 설정에 따라 이메일을 받을 수 있습니다.
How to generate a random complex unitary matrix whose columns each sum up to 1
조회 수: 51 (최근 30일)
이전 댓글 표시
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
2020년 11월 22일
편집: 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
Bruno Luong
2020년 11월 22일
편집: Bruno Luong
2020년 11월 22일
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
2020년 11월 22일
편집: David Goodmanson
2020년 11월 22일
Hi Bruno, yeah, I was feeling the pressure (in a good way) because I figured you were probably on the job.
Paul
2020년 11월 22일
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
2020년 11월 22일
편집: Bruno Luong
2020년 11월 22일
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
2020년 11월 22일
편집: Paul
2020년 11월 22일
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
2020년 11월 22일
"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
2020년 11월 23일
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
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).
Paul
2020년 11월 23일
David,
Thanks for the update. For sure your method meets the stated requirements of the OP.
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
Bruno Luong
2020년 12월 4일
편집: Bruno Luong
2020년 12월 4일
"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
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.
추가 답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Descriptive Statistics에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!오류 발생
페이지가 변경되었기 때문에 동작을 완료할 수 없습니다. 업데이트된 상태를 보려면 페이지를 다시 불러오십시오.
웹사이트 선택
번역된 콘텐츠를 보고 지역별 이벤트와 혜택을 살펴보려면 웹사이트를 선택하십시오. 현재 계신 지역에 따라 다음 웹사이트를 권장합니다:
또한 다음 목록에서 웹사이트를 선택하실 수도 있습니다.
사이트 성능 최적화 방법
최고의 사이트 성능을 위해 중국 사이트(중국어 또는 영어)를 선택하십시오. 현재 계신 지역에서는 다른 국가의 MathWorks 사이트 방문이 최적화되지 않았습니다.
미주
- América Latina (Español)
- Canada (English)
- United States (English)
유럽
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
아시아 태평양
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)