有没有更好的方法来随机生成双重随机矩阵？

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 的非零值会更有效。