title: MATLAB: Fast creation of random symmetric Matrix with fixed degree (sum of rows)
body: <p>I am searching for a method to create, in a fast way a random matrix <code>A</code> with the follwing properties:</p>
  
  <ul>
  <li><code>A = transpose(A)</code></li>
  <li><code>A(i,i) = 0</code> for all i</li>
  <li><code>A(i,j) &gt;= 0</code> for all i, j</li>
  <li><code>sum(A) =~ degree</code>; the sum of rows are randomly distributed by a distribution I want to specify (here <code>=~</code> means approximate equality).</li>
  </ul>
  
  <p>The distribution <code>degree</code> comes from a matrix <code>orig</code>, specifically <code>degree=sum(orig)</code>, thus I know that matrices with this distribution exist.</p>
  
  <p>For example: <code>orig=[0 12 7 5; 12 0 1 9; 7 1 0 3; 5 9 3 0]</code></p>
  
  <pre><code>orig =
   0    12     7     5
  12     0     1     9
   7     1     0     3
   5     9     3     0
  
   sum(orig)=[24 22 11 17];
  </code></pre>
  
  <p>Now one possible matrix <code>A=[0 11  5  8, 11  0  4  7, 5  4  0  2, 8  7  2  0]</code> is </p>
  
  <pre><code>A = 
   0    11     5     8
  11     0     4     7
   5     4     0     2
   8     7     2     0
  </code></pre>
  
  <p>with <code>sum(A)=[24 22 11 17]</code>.</p>
  
  <p>I am trying this for quite some time, but unfortunatly my two ideas didn't work:</p>
  
  <hr>
  
  <p><strong>version 1:</strong></p>
  
  <p>I switch <code>Nswitch</code> times two random elements: <code>A(k1,k3)--; A(k1,k4)++; A(k2,k3)++; A(k2,k4)--;</code> (the transposed elements aswell).</p>
  
  <p>Unfortunatly, <code>Nswitch = log(E)*E</code> (with <code>E=sum(sum(nn))</code>) in order that the Matrices are very uncorrelated. As my <code>E &gt; 5.000.000</code>, this is not feasible (in particular, as I need at least 10 of such matrices).</p>
  
  <hr>
  
  <p><strong>version 2:</strong></p>
  
  <p>I create the matrix according to the distribution from scratch. The idea is, to fill every row <code>i</code> with <code>degree(i)</code> numbers, based on the distribution of <code>degree</code>:</p>
  
  <pre><code>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;
  </code></pre>
  
  <p>However, while <code>sum(A')=degree</code>, <code>sum(A)</code> systematically deviates from <code>degree</code>, and I am not able to find the reason for that.</p>
  
  <p>Small deviations from <code>degree</code> are fine of course, but there seem to be systmatical deviations in particulat of the matrices contain in some places large numbers.</p>
  
  <hr>
  
  <p>I would be very happy if somebody could either show me a fast method for version1, or a reason for the systematic deviation of the distribution in version 2, or a method to create such matrices in a different way. Thank you!</p>
  
  <hr>
  
  <p><strong>Edit:</strong></p>
  
  <p>This is the problem in <em>matsmath</em>'s proposed solution:
  Imagine you have the matrix:</p>
  
  <pre><code>orig = 
   0    12     3     1
  12     0     1     9
   3     1     0     3
   1     9     3     0
  </code></pre>
  
  <p>with <code>r(i)=[16 22 7 13]</code>.</p>
  
  <ul>
  <li>Step 1: r(1)=16, my random integer partition is p(i)=[0 7 3 6].</li>
  <li>Step 2: Check that all p(i)&lt;=r(i), which is the case.</li>
  <li>Step 3:</li>
  </ul>
  
  <p>My random matrix starts looks like</p>
  
  <pre><code>A = 
   0    7     3     6
   7    0     .     .
   3    .     0     .
   6    .     .     0
  </code></pre>
  
  <p>with the new row sum vector rnew=[r(2)-p(2),...,r(n)-p(n)]=[15 4 7]</p>
  
  <p>Second iteration (here the problem occures):</p>
  
  <ul>
  <li>Step 1: rnew(1)=15, my random integer partition is p(i)=[0 A B]: rnew(1)=15=A+B.</li>
  <li>Step 2: Check that all p(i)&lt;=rnew(i), which gives A&lt;=4, B&lt;=7. So A+B&lt;=11, but A+B has to be 15. <strong>contradiction</strong> :-/</li>
  </ul>
  
  <hr>
  
  <p><strong>Edit2:</strong></p>
  
  <p>This is the code representing (to the best of my knowledge) the solution posted by <em>David Eisenstat</em>:</p>
  
  <pre><code>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]
  </code></pre>