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

Question

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

algorithmperformancematlabmatrix

7

我正在寻找一种快速方法来创建一个随机矩阵 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]

- Mario Krenn

sum(A) =~ degree 这句话的意思是你想让 sum(A) 成为 degree 的一个排列，还是这只是近似相等？ - David Eisenstat

抱歉造成困惑，它的意思是“近似相等”。 - Mario Krenn

矩阵是否固定为4x4？还是您正在寻求适用于所有矩阵尺寸的解决方案？ - mban

它应该适用于任意大小，我的最终矩阵将是约为5.000^2到10.000^2。因此速度在最后也是一个问题。 - Mario Krenn

谢谢tvo。这确实是一个非常相似的问题，您在那里的答案很有趣。但是我想知道如何满足所有限制条件以使您的辅助矩阵最终为非负数矩阵？对我来说，这个问题看起来与原始问题同样复杂 - 但也许我错过了一些明显的点。谢谢！ - Mario Krenn

显示剩余3条评论

2个回答

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_i，我们就可以使用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'。其逆矩阵由 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)时间内完成。我会在有机会时添加更多细节。

- David Eisenstat

谢谢，这听起来像是一个有趣而且非常快速的解决方案。不过有两个问题：我想要A（i，i）= 0，在你的方法中会有问题吗？如果我之后将对角线设置为零，那么分布可能会发生严重改变，对吧？另外，您能否向我解释一下为什么我们希望权重满足w_i*w_j？谢谢 - Mario Krenn

@NicoDean 第一个问题：没有问题，这就是为什么我建议使用梯度下降而不是只给出闭式解的起始解决方案。第二个问题：我猜测，对于满足约束条件的随机矩阵，这些将是条目的预期值。为什么是泊松分布？这有点像一遍又一遍地选择随机边。 - David Eisenstat

谢谢提供的信息。我现在已经实现了您的想法，使用上面的矩阵在我的编辑中（从“0, 12, 3, 1”开始）。解方程组给出了w_i=[2.2406 4.6334 0.8174 1.6902]。不幸的是，生成的矩阵与真实度数存在巨大的系统偏差：sum(orig)=[16, 22, 7, 13]，sum(RM)=[13.5230, 28.8452, 4.9635, 10.6683]，sum(RM')=[14.8337, 9.6171, 18.0627, 15.4865]，其中RM是使用您的算法生成的1000个随机矩阵的平均值。（如果我再运行1000次，总和非常相似，因此这是一个系统效应。:-/） - Mario Krenn

@NicoDean 应该是

a=poissrnd(w(1)*w(2));b=poissrnd(w(1)*w(3));c=poissrnd(w(1)*w(4));d=poissrnd(w(2)*w(3));e=poissrnd(w(2)*w(4));f=poissrnd(w(3)*w(4));rndmat=[0,a,b,c;a,0,d,e;b,d,0,f;c,e,f,0]

。 - David Eisenstat

@NicoDean 我很好奇你是否能在这里发布你的fsolve实现。我正在尝试做一些非常相似的事情（虽然是使用未加权的二进制图），但是在通过fsolve复制你的权重结果方面遇到了一些困难（我正在通过手动牛顿实现接近）。另外，如果您对速度优化有任何见解，也很感兴趣，这是需要重复100次以上的归一化过程，每个过程中有超过10k个排列组合（针对一系列更大的度分布）。 - yatakaka

显示剩余5条评论

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Matsmath · Accepted Answer

第一种方法（基于version2）

让你的行求和向量由矩阵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)。

也就是说，任何行总和，比如第i个，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和r3，以及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，而第一种方法根本不使用这样的知识。