Help with vectorizing Diag command
2 views (last 30 days)
Show older comments
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
0 Comments
Accepted Answer
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
0 Comments
More Answers (3)
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
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.
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;
0 Comments
See Also
Categories
Find more on Logical 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!