MATLAB：快速创建具有固定度数（行总和）的随机对称矩阵

Question

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

• A = transpose(A)
• A(i,i) = 0 为所有我
• 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 的分布，用degree(i) 数字填充每一行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)
• 第三步：

我的随机矩阵开始看起来像

A = 
 0    7     3     6
 7    0     .     .
 3    .     0     .
 6    .     .     0

使用新的行和向量 rnew=[r(2)-p(2),...,r(n)-p(n)]=[15 4 7]

第二次迭代（这里出现问题）：

• 第一步：rnew(1)=15，我的随机整数分区是p(i)=[0 A B]：rnew(1)=15=A+B。
• 第 2 步：检查所有 p(i)矛盾 :-/

---

编辑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]

Answer 1

由于大致保留度数序列就足够了，让我提出一个随机分布，其中对角线上方的每个条目都是根据Poisson distribution 选择的。我的直觉是我们想要找到权重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.

后者可以通过应用Newton's method 来解决，如下所述，从w_i = d_i / sqrt(sum_j d_j) 的初始解决方案开始。

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

---

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

4乘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'。逆由Sherman--Morrison 公式给出

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

对于牛顿步，我们需要为一些给定的y 计算J^-1*y。这可以使用上面的公式在时间O(n) 中直接完成。有机会我会补充更多细节。

Answer 2

第一种方法（基于版本2）

让你的行求和向量由矩阵 orig [r(1),r(2),...,r(n)] 给出。

Step 1. Take a random integer partition of the integer r(1) into exactly n-1 parts, say p(2), p(3), ..., p(n)
Step 2. Check if p(i)<=r(i) for all i=2...n. If not, go to Step 1.
Step 3. Fill out your random matrix first row and colum by the entries 0, p(2), ... , p(n), and consider the new row sum vector [r(2)-p(2),...,r(n)-p(n)].

对 n-1 阶矩阵重复这些步骤。

关键是，您一次随机化一行，然后将问题简化为搜索大小为小一的矩阵。

---

正如 OP 在评论中指出的那样，这种幼稚的算法失败了。原因是所讨论的矩阵在其条目上具有进一步必要条件，如下所示：

事实：

如果A 是一个具有行和[r(1), r(2), ..., r(n)] 的原矩阵，那么必然对于每个i=1..n 它都持有r(i)<=-r(i)+sum(r(j),j=1..n)。

也就是说，任何行总和，比如ith、r(i)，必然最多与其他行总和的总和一样大（不包括r(i)）。

鉴于此，修改算法是可能的。请注意，在步骤 2b 中。我们检查新的行和向量是否具有上面讨论的属性。

Step 1. Take a random integer partition of the integer r(1) into exactly n-1 parts, say p(2), p(3), ..., p(n)
Step 2a. Check if p(i)<=r(i) for all i=2...n. If not, go to Step 1.
Step 2b. Check if r(i)-p(i)<=-r(i)+p(i)+sum(r(j)-p(j),j=2..n) for all i=2..n. If not, go to Step 1.
Step 3. Fill out your random matrix first row and colum by the entries 0, p(2), ... , p(n), and consider the new row sum vector [r(2)-p(2),...,r(n)-p(n)].

---

第二种方法（基于版本1）

我不确定这种方法是否给您随机矩阵，但它肯定给您不同矩阵。

这里的想法是在本地更改 orig 矩阵的某些部分，以保持其所有属性的方式。

您应该在主对角线下方寻找一个随机的 2x2 子矩阵，其中包含严格的正条目，例如 [[a,b],[c,d]]，并通过随机值 r 到 [[a+r,b-r],[c-r,d+r]] 扰乱其内容。您也在主对角线上方进行相同的更改，以保持新矩阵对称。这里的重点是条目中的更改相互“抵消”。

当然，r 的选择方式应该是 b-r>=0 和 c-r>=0。

您也可以按照这个想法来修改更大的子矩阵。例如，您可以选择 3 个随机行坐标 r1、r2、r2 和 3 个随机列坐标 c1、c2 和 c3，然后在 orig 矩阵中进行更改9 个位置(ri,cj) 如下：您将3x3 子矩阵[[a b c],[d e f], [g h i]] 更改为[[a-r b+r c] [d+r e f-r], [g h-r i+r]]。你在转置的地方做同样的事情。同样，随机值r 的选择方式必须使a-r>=0 和f-r>=0 和h-r>=0。此外，c1 和 r1，c3 和 r3 必须是不同的，因为您不能更改矩阵 orig 的主对角线中的 0 项。

你可以一遍又一遍地重复这样的事情，比如100 次，直到你找到看起来随机的东西。请注意，这个想法使用了您已经了解解决方案的事实，这是矩阵orig，而第一种方法根本不使用此类知识。