Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Block diagonalize rows of a J-by-2 matrix

Tags:

matlab

Given a J-by-2 matrix, for example

A = [1 2 ; 3 4 ; 5 6]

I want to block diagonalize it. That is, I want:

B = [1 2 0 0 0 0 ; 0 0 3 4 0 0 ; 0 0 0 0 5 6].

One command that does this is:

blkdiag(A(1,:),A(2,:),A(3,:))

This is will be slow and tedious if J is large. Is there a built-in Matlab function that does this?

like image 657
Matthew Thirkettle Avatar asked Dec 15 '22 10:12

Matthew Thirkettle


1 Answers

Here's one hacky solution for J x 2 array case using linear indexing -

%// Get number of rows
N = size(A,1);                   

%// Get linear indices of the first column elements positions in output array
idx = 1:2*N+1:(N-1)*(2*N+1)+1;   

%// Setup output array
out = zeros(N,N*2);

%// Put first and second column elements into idx and idx+N positions
out([idx(:) idx(:)+N]) = A

There is just one function call (neglecting size as it must be minimal) overhead of zeros and even that can be removed with this undocumented zeros initialization trick -

out(N,N*2) = 0; %// Instead of out = zeros(N,N*2);

Sample run -

A =
     1     2
     3     4
     5     6
     7     8
out =
     1     2     0     0     0     0     0     0
     0     0     3     4     0     0     0     0
     0     0     0     0     5     6     0     0
     0     0     0     0     0     0     7     8

Here's a benchmarking for the posted solutions thus far.

Benchmarking Code

%//Set up some random data
J = 7000;   A = rand(J,2);
%// Warm up tic/toc
for k = 1:100000
    tic(); elapsed = toc();
end   
disp('---------------------------------- With @mikkola solution')
tic
temp = mat2cell(A, ones(J,1), 2);
B = blkdiag(temp{:});
toc, clear B temp
disp('---------------------------------- With @Jeff Irwin solution')
tic
m = size(A, 1);
n = size(A, 2);
B = zeros(m, m * n);
for k = 1: n
    B(:, k: n: m * n) = diag(A(:, k));
end
toc, clear B k m n
disp('---------------------------------- With Hacky1 solution')
tic
N = size(A,1);                   
idx = 1:2*N+1:(N-1)*(2*N+1)+1;   
out = zeros(N,N*2);
out([idx(:) idx(:)+N]) = A;
toc, clear out idx N
disp('---------------------------------- With Hacky2 solution')
tic
N = size(A,1);                   
idx = 1:2*N+1:(N-1)*(2*N+1)+1;
out(N,N*2) = 0;
out([idx(:) idx(:)+N]) = A;
toc, clear out idx N

Runtimes

---------------------------------- With @mikkola solution
Elapsed time is 0.546584 seconds.
---------------------------------- With @Jeff Irwin solution
Elapsed time is 1.330666 seconds.
---------------------------------- With Hacky1 solution
Elapsed time is 0.455735 seconds.
---------------------------------- With Hacky2 solution
Elapsed time is 0.364227 seconds.
like image 198
Divakar Avatar answered Dec 28 '22 04:12

Divakar