MATLAB Answers

Efficient algorithm for a duplication matrix

48 views (last 30 days)
Youngkyu Kim
Youngkyu Kim on 27 Jul 2019
Commented: Youngkyu Kim on 30 Jul 2019
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(:))

  2 Comments

Walter Roberson
Walter Roberson on 27 Jul 2019
What are vec and vech in this context?

Sign in to comment.

Accepted Answer

Jan
Jan on 28 Jul 2019
Edited: Jan on 29 Jul 2019
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
r = r + n - i + 1;
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! :-)

  2 Comments

Jan
Jan on 30 Jul 2019
@Youngkyu Kim: The runtime behavior of the original code is quadratic and the extrapolated run-time for n=1000 is about 14'000 seconds. My last suggestions needs 0.12 seconds.
I took me some time to figure this out. Don't you think, that the problem is solved? A short reaction would be fine, or at least accepting the answer as a working solution. I do not need any reputation points anymore, but I think, such a success is worth to be honored.
Youngkyu Kim
Youngkyu Kim on 30 Jul 2019
Thanks for your suggesteion.
I modified my code and it worked well enough.
Your suggestion is slightly faster than mine.
function D = duplication(n)
m=1/2*n*(n+1);
nsq=n^2;
Lind=tril(true(n));
Lind=find(Lind(:));
Lind=Lind(:);
Uind=rem(Lind-1,n)*n+ceil(Lind/n);
i=(1:m)';
a=(i-1)*nsq+Lind;
b=(i-1)*nsq+Uind;
c=union(a,b);
[I,J]=ind2sub([nsq,m],c);
D=sparse(I,J,1,nsq,m);
end

Sign in to comment.

More Answers (0)

Sign in to answer this question.