Asked by Youngkyu Kim
on 27 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(:))

Answer by Jan
on 28 Jul 2019

Edited by Jan
on 29 Jul 2019

Accepted 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

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! :-)

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
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.

Opportunities for recent engineering grads.

Apply Today
## 2 Comments

## Walter Roberson (view profile)

## Direct link to this comment

https://se.mathworks.com/matlabcentral/answers/473737-efficient-algorithm-for-a-duplication-matrix#comment_729028

## Stephan (view profile)

## Direct link to this comment

https://se.mathworks.com/matlabcentral/answers/473737-efficient-algorithm-for-a-duplication-matrix#comment_729030

Sign in to comment.