MATLAB:快速创建具有固定度(行总和)的随机对称矩阵

Question

我正在寻找一种方法来快速创建A具有以下属性的随机矩阵:

• A = transpose(A)
• A(i,i) = 0 对于所有我
• A(i,j) >= 0 对于所有我,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:

我从头开始根据分布创建矩阵.这个想法是,根据以下分布填充每一行i的degree(i)数字degree:

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 AB]:rnew(1)= 15 = A + B.
• 步骤2:检查所有p(i)<= rnew(i),得到A <= 4,B <= 7.所以A + B <= 11,但A + B必须是15. 矛盾: - /

---

EDIT2:

这是代表(据我所知)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]

Answer 1

由于近似保留度数序列就足够了，所以让我提出一个随机分布，其中对角线上方的每个条目都是根据泊松分布选择的。我的直觉是，我们想要找到权重w_i，使得i,j的条目i != j具有均值w_i*w_j（所有对角线条目均为零）。这给了我们一个非线性方程组：

for all i, (sum_{j != i} w_i*w_j) = d_i,

其中d_i是 的度i。等价地，

for all i, w_i * (sum_j w_j) - w_i^2 = d_i.

后者可以通过应用如下所述的牛顿法w_i = d_i / sqrt(sum_j d_j)从 的起始解来求解。

一旦我们有了w_is，我们就可以使用重复采样poissrnd来一次生成多个泊松分布的样本。

（如果我有时间，我会尝试在 numpy 中实现这个。）

4 x 4 问题方程组的雅可比矩阵为

(w_2 + w_3 + w_4) w_1               w_1               w_1
w_2               (w_1 + w_3 + w_4) w_2               w_2
w_3               w_3               (w_1 + w_2 + w_4) w_3
w_4               w_4               w_4               (w_1 + w_2 + w_3).

一般来说，令A为对角矩阵，其中A_{i,i} = sum_j w_j - 2*w_i. 让u = [w_1, ..., w_n]'和v = [1, ..., 1]'. 雅可比行列式可以写成J = A + u*v'。倒数由谢尔曼-莫里森公式给出

                              A^-1*u*v'*A^-1
J^-1 = (A + u*v')^-1 = A^-1 - -------------- .
                              1 + v'*A^-1*u

对于牛顿步，我们需要计算J^-1*y一些给定的y。O(n)这可以使用上面的方程直接及时完成。当我有机会时，我会添加更多细节。