双随机软件java,有没有更好的方法来随机生成双随机矩阵？

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

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

% 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

解决方案

How about changing lines 14 and 15 to the following lines:

l = ( [ pidx ; 1:n ] - 1 ) * [1;n] + 1; % convert pairs (pidx,1:n) to linear indices

M(l) = M(l) + c(k);

since P is very sparse, maybe it would be more efficient to increment only the non-zeros of P.