## Is there a better way to randomly generate a Doubly Stochastic Matrix?

on 17 Nov 2012

### Matt Fig (view profile)

Here is a profile of the call n = 2*10^3; M = DStochMat02(n,ones(n)./n);

More specifically, can the hot-spot, statement 14, be reduced?

```   time   calls  line
1 function M = DStochMat02(n,c)
2 % Generate a random doubly stochastic matrix using
3 % Theorem (Birkhoff [1946], von Neumann [1953])
4 % Any doubly stochastic matrix M can be written as a     convex combination
5 % of permutation matrices P1,...,Pk (i.e. M = c1*P1+...+ ck*Pk
6 % for nonnegative c1,...,ck with c1+...+ck = 1).
7 % Complexity: O(n^2)
8 % USE: M = DStochMat02(4,[1/2 1/8 1/8 1/4])
9 % Derek O'Connor, Oct 2006, Nov 2012. derekroconnor@eircom.net```
```    0.02       1   10 M = zeros(n,n);
< 0.01       1   11 I = eye(n);
< 0.01       1   12 for k = 1:n
1.64   2000   13     pidx = GRPdur(n);  % Random Permutation
107.72   2000   14     P = I(pidx,:);     % Random P matrix
41.09   2000   15     M = M + c(k)*P;
< 0.01    2000   16 end```
```%--------------------------------------------------------------
function p = GRPdur(n)
% -------------------------------------------------------------
% Generate a random permutation p(1:n) using Durstenfeld's
% O(n) Shuffle Algorithm, CACM, 1964.
% See Knuth, Section 3.4.2, TAOCP, Vol 2, 3rd Ed.
% Complexity: O(n)
% USE: p = GRPdur(10^7);
% Derek O'Connor, 8 Dec 2010.  derekroconnor@eircom.net
% -------------------------------------------------------------
```
```      p = 1:n;                  % Start with Identity permutation
for k = n:-1:2
r = 1+floor(rand*k);      % random integer between 1 and k
t    = p(k);
p(k) = p(r);               % Swap(p(r),p(k)).
p(r) = t;
end
return % GRPdur```

Matt Fig

### Matt Fig (view profile)

on 17 Nov 2012

Can you post what GRPdur does? Is it akin to RANDPERM?

## Products

### Matt Fig (view profile)

on 17 Nov 2012
Edited by Matt Fig

### Matt Fig (view profile)

on 17 Nov 2012

This looks promising. I have kept some notes in the comments so you can see what I was doing...

There are three lines of code that are commented out. If you uncomment them, and comment the line after the second one, you can check that these give the same results. Of course the timing is not fair under those circumstances, but with the code commented the timings show a fantastic speed difference for large values....

```M = zeros(n,n);
I = eye(n);
```
```tic  % Original method.
for k = 1:n
pidx = GRPdur(n);  % Random Permutation
%     F{k} = pidx;     % Used to check same results below....
P = I(pidx,:);     % Random P matrix
M = M + c(k)*P;
end
toc
```
```M2 = zeros(n,n);
```
```tic  % proposed alternative
for k = 1:n
%     [~,idx] = sort(F{k});    % Used to compare same results.
[~,idx] = sort(GRPdur(n)); % Or just idx = GRPdur(n);  ??
idx = idx + (0:n-1)*n;
M2(idx) = M2(idx) + c(k);
end
toc
```
```% max(abs(M2(:)-M(:)))  % Used to compare same results.
```

Matt Fig

### Matt Fig (view profile)

on 18 Nov 2012

@Matt, Amazing, obscure, and very fast (x70), and it seems to give the correct results. I need to do more testing.

Now, can you explain, mathematically, what you are doing? You have transformed the problem somehow. Can you explain?

```function times = TestStochMatFig(nvals);
%
nalgs = 2;
times = zeros(length(nvals),nalgs);
kn = 0;
for n = nvals
kn = kn+1;
M = zeros(n,n);
I = eye(n);
c = rand(n,1);
c = c./sum(c);
algno = 1;
M1 = M;
tic;  % Original method. O'Connor
for k = 1:n
pidx = GRPdur(n);  % Random Permutation
P = I(pidx,:);     % Random P matrix
M1 = M1 + c(k)*P;
end
times(kn,algno) = toc;
algno = 2;
M2 = M;
tic;   % proposed alternative. Matt Fig
for k = 1:n
[~,idx] = sort(GRPdur(n)); % Or just idx = GRPdur(n);  ??
idx = idx + (0:n-1)*n;
M2(idx) = M2(idx) + c(k);
end
times(kn,algno) = toc;
max(abs(M2(:)-M1(:)))  % Used to compare same results.
```
```end
```
Matt Fig

### Matt Fig (view profile)

on 18 Nov 2012

Hi Derek, I am glad you found the results useful. What I did was basically take the operations from matrix operations to vector operations. Since you are multiplying a scalar (c(k)) by an identity matrix, there is no need to keep all the zeros. I simply take the single vector of values (there are n c(k) values) and index that directly into only those elements of M that are actually updating - thus avoiding the unnecessary addition of

```K = n^2 - n
```

zeros each time through the loop. There is no need to add K zeros to M each time through! In addition, my approach avoids multiplying K+n numbers by c(k) each time through the loop, when we don't even need to do a single such multiplication. When n is large these operations takes a significant amount of time and memory acrobatics....

Also, note that I only kept the call to SORT in there in order to make sure that the codes produced the exact same result (when the appropriate lines are uncommented of course). After hundreds of runs, my code does produce the exact same result. Thus I think you could avoid the call to SORT all together unless you have some statistical connection to GRPdur you are trying to preserve.

### Derek O'Connor (view profile)

on 23 Nov 2012

@Matt, You are a true blue Matlabber. With just matrix indexing you have written a function which is many times faster than mine. See below

```    % -------------------------------------------------------------
function table = TestStochMat02Fig(nvals);
% -------------------------------------------------------------
% Testing two functions for generating random doubly stochastic
% matrices
% USE: times = TestStochMat02Fig([10^3:10^3:6*10^3]);
% Derek O'Connor, Matt Fig, 20 Nov 2012
% -------------------------------------------------------------
nalgs = 2;
times = zeros(length(nvals),nalgs);
kn = 0;
for n = nvals
kn = kn+1;
M = zeros(n,n);
I = eye(n);
c = rand(n,1);
c = c./sum(c);
%---------------- Function by Derek O'Connor -------------
algno = 1;
M1 = M;
tic;
for k = 1:n
p = GRPdur(n);                  % Random Permutation p
P = I(p,:);                     % Random P-matrix P
M1 = M1 + c(k)*P;
end
times(kn,algno) = toc;
%---------------- Function by Matt Fig--------------------
algno = 2;
M2 = M;
tic;
for k = 1:n
p = GRPdur(n);
p = p + (0:n-1)*n;
M2(p) = M2(p) + c(k);
end
times(kn,algno) = toc;```
```        max(abs(M2(:)-M1(:)))  % Used to compare same results.
end % n = nvals
table = [nvals' times(:,1) times(:,2) times(:,1)./times(:,2)];   ```
```% ----------- End: table = TestStochMat02Fig(nvals); -------------
```
```>table = TestStochMat02Fig([10^3:10^3:6*10^3]);
```
```                n         Derek        Matt         D/M
1000       6.6936      0.10741       62.317
2000       60.822      0.51116       118.99
3000       209.14       1.1846       176.54
4000       453.56       1.8474       245.51
5000       897.59       3.7354       240.29
6000       1490.2       4.8622       306.49```

Be warned! There is a trickier problem to come.