Help with vectorizing Diag command

5 views (last 30 days)
Brendan
Brendan on 7 Jul 2012
I have an MxN matrix, which I would like to turn into an MxNxN matrix, where the elements in the second and third dimensions are a diagonal matrix. It can be done in a for loop by
for ii=1:M
newMat(ii,:,:)=diag(oldMat(ii,:))
end
Any ideas on how to vectorize this? Thanks

Accepted Answer

Jan
Jan on 7 Jul 2012
Edited: Jan on 7 Jul 2012
Pre-allocate:
clear
M = 200;
N = 300;
oldMat = rand(M, N);
tic; % ORIGINAL
for ii=1:M
newMat(ii,:,:) = diag(oldMat(ii,:));
end
toc % 16.71 seconds (Matlab 2009a/64, Win7, Core2Duo)
tic; % REPMAT
[m, n] = size(oldMat);
idx = permute(repmat(logical(eye(n)), [1, 1, m]), [3, 2, 1]);
nM = idx+0;
nM(idx) = oldMat;
toc % 1.10 sec
clear newMat
tic; % WITH PRE-ALLOCATION
newMat = zeros(M, N, N); % Pre-allocate!!!
for ii=1:M
newMat(ii,:,:) = diag(oldMat(ii,:));
end
toc % 0.68 seconds
clear newMat
tic; % BSXFUN
newMat = zeros(M,N,N);
newMat(bsxfun(@plus,(1:M)', (0:N-1)*(M*N+M))) = oldMat;
toc % 0.076 seconds

More Answers (3)

Sven
Sven on 7 Jul 2012
Edited: Sven on 7 Jul 2012
Hi Brendan, is this the answer you're looking for?
M = 2, N = 3
oldMat = reshape(1:M*N,M,N)
newMat = zeros(M,N,N)
newMat(bsxfun(@plus,(1:M)',(0:N-1)*M*(N+1))) = oldMat
It uses the bsxfun command to build a list of indices equivalent to the non-zero elements of your newMat, and just populates these elements directly from oldMat.
I haven't tested, but I suspect it will be faster than the loop. Hope it's helpful.
[Edit]
Just tested for M = 200, N = 300:
Elapsed time is 8.962171 seconds. <- loop method
Elapsed time is 1.915089 seconds. <- bsxfun method
  2 Comments
Jan
Jan on 7 Jul 2012
In my measurements BSXFUN is much faster.
Sven
Sven on 7 Jul 2012
Oh, I wasn't timing just once... I gave it a few loops which is why my numbers are over 1 second :) Either that or my pc runs like treacle.
And for fairness to the loop method, my comparison used preallocation.

Sign in to comment.


Andrei Bobrov
Andrei Bobrov on 7 Jul 2012
[m,n] = size(oldMat);
idx = permute(repmat(logical(eye(n)),[1,1,m]),[3 2 1]);
nM = idx+0;
nM(idx)=oldMat;

Brendan
Brendan on 11 Jul 2012
Thanks a lot!

Community Treasure Hunt

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

Start Hunting!