有沒有更好的方法來隨機生成雙隨機矩陣？

Question

這裏是調用的輪廓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 

Answer 1

如何改變線路14和15以下行：

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

因爲P是很稀疏，也許這將是更有效地增加只有P的非零。