Efficient algorithm for a duplication matrix

16 views (last 30 days)
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(:))

Accepted Answer

Jan
Jan on 28 Jul 2019
Edited: Jan on 4 Aug 2021
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);
cn = cumsum(n:-1:2); % [EDITED, 2021-08-04], 10% faster
for i = 1:n
% v(r:r + i - 2) = i - n + cumsum(n - (0:i-2));
v(r:r + i - 2) = i - n + cn(1:i - 1); % [EDITED, 2021-08-04]
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! :-)
  7 Comments
Jan
Jan on 5 Aug 2021
The list C of my code needs 6GB of RAM, if all 1000 matrices are created. Writing it as -v7.3 MAT file creates an 1GB file. Reading it takes 15 seconds on a HDD instead of 35 seconds of creating all matrices dynamically.

Sign in to comment.

More Answers (0)

Categories

Find more on Loops and Conditional Statements in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!