我正在寻找一种快速方法来创建一个随机矩阵 A
,并满足以下属性:
A = transpose(A)
A(i,i) = 0
,其中 i 是任意整数A(i,j) >= 0
,其中 i 和 j 是任意整数sum(A) =~ degree
;每行的和按照我指定的分布进行随机分布(这里的=~
表示近似相等)。
分布 degree
来自于矩阵 orig
,具体来说是 degree=sum(orig)
,因此我知道具有这种分布的矩阵是存在的。
例如:orig=[0 12 7 5; 12 0 1 9; 7 1 0 3; 5 9 3 0]
orig =
0 12 7 5
12 0 1 9
7 1 0 3
5 9 3 0
sum(orig)=[24 22 11 17];
现在可能有一个矩阵A=[0 11 5 8, 11 0 4 7, 5 4 0 2, 8 7 2 0]
A =
0 11 5 8
11 0 4 7
5 4 0 2
8 7 2 0
使用sum(A)=[24 22 11 17]
。
我已经尝试了很长时间,但不幸的是我的两个想法都没有奏效:
版本1:
我交换Nswitch
次两个随机元素:A(k1,k3)-;A(k1,k4)+;A(k2,k3)+;A(k2,k4)-;
(同样交换转置的元素)。
不幸的是,为了使矩阵非常不相关,Nswitch = log(E)*E
(其中E=sum(sum(nn))
)。因为我的E > 5,000,000
,所以这是不可行的(特别是因为我需要至少10个这样的矩阵)。
版本2:
我根据分布从头开始创建矩阵。我的想法是,根据degree
的分布填充每一行i
中的degree(i)
个数字:
nn=orig;
nnR=zeros(size(nn));
for i=1:length(nn)
degree=sum(nn);
howmany=degree(i);
degree(i)=0;
full=rld_cumsum(degree,1:length(degree));
rr=randi(length(full),[1,howmany]);
ff=full(rr);
xx=i*ones([1,length(ff)]);
nnR = nnR + accumarray([xx(:),ff(:)],1,size(nnR));
end
A=nnR;
然而,虽然
sum(A')=degree
,但sum(A)
与degree
有系统偏差,我无法找到原因。当然,对于
degree
的小偏差是可以接受的,但在某些包含大数值的矩阵中存在系统性偏差。
如果有人能向我展示快速解决方案版本1,或者给出版本2分布系统偏差的原因,或者提供一种不同的创建这种矩阵的方法,我将非常感激。谢谢!
编辑:
matsmath提出的解决方案存在以下问题:
想象一下你有一个矩阵:
orig =
0 12 3 1
12 0 1 9
3 1 0 3
1 9 3 0
使用r(i)=[16 22 7 13]
。
- 步骤1: r(1)=16,我的随机整数分区是p(i)=[0 7 3 6]。
- 步骤2: 检查所有的p(i)<=r(i),这是成立的。
- 步骤3:
我的随机矩阵开始看起来像
A =
0 7 3 6
7 0 . .
3 . 0 .
6 . . 0
使用新的行和向量 rnew=[r(2)-p(2),...,r(n)-p(n)]=[15 4 7]
第二次迭代(问题出现在这里):
- 步骤1:rnew(1)=15,我的随机整数分区是 p(i)=[0 A B]:rnew(1)=15=A+B。
- 步骤2:检查所有 p(i)<=rnew(i),这给出 A<=4,B<=7。所以 A+B<=11,但 A+B 必须为15。矛盾 :-/
编辑2:
这是代表(据我所知)David Eisenstat 发布的解决方案的代码:
orig=[0 12 3 1; 12 0 1 9; 3 1 0 3; 1 9 3 0];
w=[2.2406 4.6334 0.8174 1.6902];
xfull=zeros(4);
for ii=1:1000
rndmat=[poissrnd(w(1),1,4); poissrnd(w(2),1,4); poissrnd(w(3),1,4); poissrnd(w(4),1,4)];
kkk=rndmat.*(ones(4)-eye(4)); % remove diagonal
hhh=sum(sum(orig))/sum(sum(kkk))*kkk; % normalisation
xfull=xfull+hhh;
end
xf=xfull/ii;
disp(sum(orig)); % gives [16 22 7 13]
disp(sum(xf)); % gives [14.8337 9.6171 18.0627 15.4865] (obvious systematic problem)
disp(sum(xf')) % gives [13.5230 28.8452 4.9635 10.6683] (which is also systematically different from [16, 22, 7, 13]
sum(A) =~ degree
这句话的意思是你想让sum(A)
成为degree
的一个排列,还是这只是近似相等? - David Eisenstat5.000^2
到10.000^2
。因此速度在最后也是一个问题。 - Mario Krenn