MATLAB: Efficient algorithm for a duplication matrix

linalgsymmetric matrixvectorization

Can anybody help me to design a Matlab code function that creates a duplication matrix D?
Thanks in advnace.
My codes is very slow…
Any ideas to speed it up?
n=1000;
% Duplication matrix: vec(P)=Dvech(P)
tic
m=1/2*n*(n+1);
nsq=n^2;
DT=sparse(m,nsq);
for j=1:n
    for i=j:n
        ijth=(j-1)*n+i;
        jith=(i-1)*n+j;
        vecTij=sparse(ijth,1,1,nsq,1);
        vecTij(jith,1)=1;
        k=(j-1)*n+i-1/2*j*(j-1);
        uij=sparse(k,1,1,m,1);
        DT=DT+uij*vecTij';
    end
end
D=DT';
toc
% test duplication matrix
C=rand(n,n);
P=1/2*(C+C');
vechP=nonzeros(tril(P));
vecP=P(:);
err_D=vecP-D*vechP;
max(err_D(:))
min(err_D(:))

Best Answer

For n=300 this needs 1.3 sec instead of 27.5 sec:
tic
m   = n * (n + 1) / 2;
nsq = n^2;
D   = spalloc(nsq, m, nsq);
row = 1;
a   = 1;
for i = 1:n
   b = i;
   for j = 0:i-2
      D(row + j, b) = 1;
      b = b + n - j - 1;
   end
   row = row + i - 1;
   
   for j = 0:n-i
      D(row + j, a + j) = 1;
   end
   row = row + n - i + 1;
   a   = a + n - i + 1;
end
toc
But it is much faster to create the index vector at first instead of accessing the sparse matrix repeatedly:
tic
m   = n * (n + 1) / 2;
nsq = n^2;
r   = 1;
a   = 1;
v   = zeros(1, nsq);
for i = 1:n
   b = i;
   for j = 0:i-2
      v(r) = b;
      b    = b + n - j - 1;
      r    = r + 1;
   end
   
   for j = 0:n-i
     v(r) = a + j;
     r    = r + 1;
   end
   % BUGFIX: Omit "r = r + n - i + 1;"  Thanks Trisha Phippard
   a = a + n - i + 1;
end
D2 = sparse(1:nsq, v, 1, nsq, m);
toc
Now I get 0.013 sec for n=300. Finally vectorize the 2 inner loops:
tic
m   = n * (n + 1) / 2;
nsq = n^2;
r   = 1;
a   = 1;
v   = zeros(1, nsq);
for i = 1:n
   v(r:r + i - 2) = i - n + cumsum(n - (0:i-2));
   r = r + i - 1;
   
   v(r:r + n - i) = a:a + n - i;
   r = r + n - i + 1;
   a = a + n - i + 1;
end
D2 = sparse(1:nsq, v, 1, nsq, m);
toc
0.011 sec. A speedup of factor 2500 for n=300. And 0.12 sec for n=1000. Nice! :-)
Related Question