7

I am searching for a method to create, in a fast way a random matrix A with the follwing properties:

  • A = transpose(A)
  • A(i,i) = 0 for all i
  • A(i,j) >= 0 for all i, j
  • sum(A) =~ degree; the sum of rows are randomly distributed by a distribution I want to specify (here =~ means approximate equality).

The distribution degree comes from a matrix orig, specifically degree=sum(orig), thus I know that matrices with this distribution exist.

For example: 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];

Now one possible matrix A=[0 11 5 8, 11 0 4 7, 5 4 0 2, 8 7 2 0] is 

A = 0 11 5 8 11 0 4 7 5 4 0 2 8 7 2 0

with sum(A)=[24 22 11 17].

I am trying this for quite some time, but unfortunatly my two ideas didn't work:

version 1:

I switch Nswitch times two random elements: A(k1,k3)--; A(k1,k4)++; A(k2,k3)++; A(k2,k4)--; (the transposed elements aswell).

Unfortunatly, Nswitch = log(E)*E (with E=sum(sum(nn))) in order that the Matrices are very uncorrelated. As my E > 5.000.000, this is not feasible (in particular, as I need at least 10 of such matrices).

version 2:

I create the matrix according to the distribution from scratch. The idea is, to fill every row i with degree(i) numbers, based on the distribution of 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;

However, while sum(A')=degree, sum(A) systematically deviates from degree, and I am not able to find the reason for that.

Small deviations from degree are fine of course, but there seem to be systmatical deviations in particulat of the matrices contain in some places large numbers.

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!

Edit:

This is the problem in matsmath's proposed solution: Imagine you have the matrix:

orig = 0 12 3 1 12 0 1 9 3 1 0 3 1 9 3 0

with r(i)=[16 22 7 13].

  • Step 1: r(1)=16, my random integer partition is p(i)=[0 7 3 6].
  • Step 2: Check that all p(i)<=r(i), which is the case.
  • Step 3:

My random matrix starts looks like

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

with the new row sum vector rnew=[r(2)-p(2),...,r(n)-p(n)]=[15 4 7]

Second iteration (here the problem occures):

  • Step 1: rnew(1)=15, my random integer partition is p(i)=[0 A B]: rnew(1)=15=A+B.
  • Step 2: Check that all p(i)<=rnew(i), which gives A<=4, B<=7. So A+B<=11, but A+B has to be 15. contradiction :-/

Edit2:

This is the code representing (to the best of my knowledge) the solution posted by 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]
Improve this question
8
  • Does sum(A) =~ degree mean that you want sum(A) to be a permutation of degree, or is this approximate equality? Commented Apr 28, 2016 at 2:42
  • Sorry for the confusion, it means "approximate equality". Commented Apr 28, 2016 at 9:50
  • Is the matrix fixed at 4x4? Or are you asking for a solution for all matrix sizes? Commented Apr 28, 2016 at 14:19
  • It should work for arbitrary sizes, my final matrix will be ~5.000^2 to 10.000^2. So speed is also an issue in the end. Commented Apr 28, 2016 at 14:24
  • As a suggestion to your version 1 solution, would randomly choosing the number to be added/subtracted to two random elements help break the correlation quicker? In other words, in your example you add/subtract 1 from each element - if you select another number at random (but is less than or equal to the total sum of the row) to add/subtract, would it require less iterations to achieve uncorrelated matrices? Commented Apr 28, 2016 at 14:53

2 Answers 2

Reset to default
1
+50

Since it's enough to approximately preserve the degree sequence, let me propose a random distribution where each entry above the diagonal is chosen according to a Poisson distribution. My intuition is that we want to find weights w_i such that the i,j entry for i != j has mean w_i*w_j (all of the diagonal entries are zero). This gives us a nonlinear system of equations:

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

where d_i is the degree of i. Equivalently,

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

The latter can be solved by applying Newton's method as described below from a starting solution of w_i = d_i / sqrt(sum_j d_j).

Once we have the w_is, we can sample repeatedly using poissrnd to generate samples of multiple Poisson distributions at once.

(If I have time, I'll try implementing this in numpy.)

The Jacobian matrix of the equation system for a 4 by 4 problem is

(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).

In general, let A be a diagonal matrix where A_{i,i} = sum_j w_j - 2*w_i. Let u = [w_1, ..., w_n]' and v = [1, ..., 1]'. The Jacobian can be written J = A + u*v'. The inverse is given by the Sherman--Morrison formula

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

For the Newton step, we need to compute J^-1*y for some given y. This can be done straightforwardly in time O(n) using the above equation. I'll add more detail when I get the chance.

Improve this answer
Sign up to request clarification or add additional context in comments.

10 Comments

Thank you, that sounds like an interesting and very fast solution. However, there are two questions: I want A(i,i)=0, is that a problem in your approach? If I set the diagonals to zero afterwards, that might in change the distribution severly, right? And secondly, could you please explain to me a bit more why we want the weights to fulfill w_i*w_j? Thanks
@NicoDean First question: no problem there, that's why I have to suggest gradient descent as opposed to just giving that starting solution as the closed form solution. Second question: I conjecture that, for a random matrix meeting the constraints, those would be the expected values of the entries. Why a Poisson? It's sort of like picking a random edge over and over.
Thanks for the infos. I implemented your idea now, using the matrix above in my Edit (starting with 0, 12, 3, 1). Solving the equation system gives me w_i=[2.2406 4.6334 0.8174 1.6902]. Unfortunatly, the produced matrices have a huge systematic deviation of the real degree: 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], where RM is the average of 1000 random matrices produced with your algorithm. (if i run it another 1000 times, the sums are very similar thus its a systematic effect. :-/
@NicoDean You're doing something wrong, then. Post your code.
@NicoDean Should be 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].
|
1

First approach (based on version2)

Let your row sum vector given by the matrix 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)].

Repeat these steps with a matrix of order n-1.

The point is, that you randomize one row at a time, and reduce the problem to searching for a matrix of size one less.

As pointed out by OP in the comment, this naive algorithm fails. The reason is that the matrices in question have a further necessary condition on their entries as follows:

FACT:

If A is an orig matrix with row sums [r(1), r(2), ..., r(n)] then necessarily for every i=1..n it holds that r(i)<=-r(i)+sum(r(j),j=1..n).

That is, any row sum, say the ith, r(i), is necessarily at most as big as the sum of the other row sums (not including r(i)).

In light of this, a revised algorithm is possible. Note that in Step 2b. we check if the new row sum vector has the property discussed above.

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)].

Second approach (based on version1)

I am not sure if this approach gives you random matrices, but it certainly gives you different matrices.

The idea here is to change some parts of your orig matrix locally, in a way which maintains all of its properties.

You should look for a random 2x2 submatrix below the main diagonal which contains strictly positive entries, like [[a,b],[c,d]] and perturbe its contents by a random value r to [[a+r,b-r],[c-r,d+r]]. You make the same change above the main diagonal too, to keep your new matrix symmetric. Here the point is that the changes within the entries "cancel" each other out.

Of course, r should be chosen in a way such that b-r>=0 and c-r>=0.

You can pursue this idea to modify larger submatrices too. For example, you might choose 3 random row coordinates r1, r2, r2 and 3 random column coordinates c1, c2, and c3 and then make changes in your orig matrix at the 9 positions (ri,cj) as follows: you change your 3x3 submatrix [[a b c],[d e f], [g h i]] to [[a-r b+r c] [d+r e f-r], [g h-r i+r]]. You do the same at the transposed places. Again, the random value r must be chosen in a way so that a-r>=0 and f-r>=0 and h-r>=0. Moreover, c1 and r1, and c3 and r3 must be distinct as you can't change the 0 entries in the main diagonal of the matrix orig.

You can repeat such things over and over again, say 100 times, until you find something which looks random. Note that this idea uses the fact that you have existing knowledge of a solution, this is the matrix orig, while the first approach does not use such knowledge at all.

Improve this answer

3 Comments

Thank you, I see that your condition "check if p(i)<=r(i)" prevents the problem of the wrong distribution in my version2. That should be very easy to adapt, i hope it's also efficient. Thanks!
I implemented the idea now, but unfortunatly it does not work. I updated by post to show a case that leads to a contradiction. Do you have a fix for the problem or a different idea? thanks
Thanks for your update. I will try now your adapted criteria. About the second idea, it is very similar my version1. Unfortunatly, literature says you need roughly log(E)*E changes to get random, uncorrelated matrices, where E are number of edges. And I have E=5.000.000, and I need at least 10 such matrices (better would be 100) - so it's unfeasible for the hardware I have.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.