Alternative to blkdiag and mat2cell functions

Hi all,

I currently have a code which takes some matrix A and stacks each row of the aforementioned matrix in to a block diagonal matrx C such that:

The code for that is:

Nc = 10 ;   % Number of rows of A
Nb = 4 ;    % Number of columns of A

A = randn(Nc, Nb) ;
B = mat2cell(A, ones(Nc, 1), Nb) ;
C = blkdiag(B{:}) ;

Generating matrix C however is extremely slow as I run a similar piece of code as shown above a couple of thousands of times. From timing the entire thing, MATLAB has pointed out that the mat2cell and blkdiag functions are by far the slowest ones, so I thought I should replace them in the following way:

Nb = 4 ;
Nc = 10 ;

A = randn(Nc, Nb) ;
C = zeros(Nc, Nb*Nc) ; % preallocate C
ridx = repmat(1:Nc, 1, Nb)' ; % row indices of where each element of A should be in C
cidx = reshape(reshape(1:Nb*Nc, Nb, Nc)', [], 1) ; % column indices of where each element of A should be in C
C(ridx, cidx) = A ;

In short, the above code simply spots the indices of where all the elements of matrix A should be allocated on matrix C, thus cidx is the row indices and ridx is the column indices.

Even though these indices are correct, this code returns an error as the matrix C(ridx, cidx) is a [Nb*Nc x Nb*Nc] matrix (or 40x40 in this case), and A is a [Nc x Nb] matrix. I did not expect that since ridx and cidx are column vectors representing indices, so I would simply expect that all the elements of A would transfer to C at the specified indices.

So how can I make the above code run correctly?

Thanks for your help in advance.

[Bruno Luong]

You clearly didn't run the code

Nc = 10 ;   % Number of rows of A
Nb = 4 ;    % Number of columns of A

A = randn(Nb, Nc) ;
B = mat2cell(A, onces(Nc, 1), Nb) ;
C = blkdiag(B{:}) ;

[KostasK]

Apologies, made a mess out of copy pasting where I changed the name of the variables from my code to the code that I posted here. I edited this accordingly.

Answer 1

[Walter Roberson]

You could consider using sparse() to put the matrix together. It is designed to take vectors of row coordinates and corresponding column coordinates and corresponding values, and assemble them into a matrix. You could then full() if you need to.

The difficulty you are encountering is that when you index on multiple non-scalar indices, then MATLAB does not use "corresponding" elements. Instead, MATLAB uses all possible combinations of indices.

M([1 2], [3 4]) = 5

would assign to M(1,3), M(2,3), M(1,4), M(2,4)

MATLAB does not have a syntax for "corresponding" index assignment. Instead what you need to do is use a method of converting the corresponding indices into linear indices and assign to there. The function designed to work convert corresponding indices into linear indices is sub2ind().

There are occasions on which it pays to skip using sub2ind() and instead use the mathematical transformations involved -- especially if you are doing several cases of assigning to sub-arrays of consistent sizes and spatial relationship. For example, if you have the linear indices IDX for a block at (R,C) and want to address something further down the diagonal, then

new_IDX = IDX + CD*number_of_rows + RD

where CD is difference in columns and RD is difference in rows .

[KostasK]

Thanks for that. I see then that linear indexing would be the way to go rather than subscripts

Answer 2

[Bruno Luong]

What about about a simple for-loop

[Nb,Nc] = size(A);
C = zeros(Nb,Nb*Nc);
for r=1:Nb
    C(r,(r-1)*Nc+(1:Nc)) = A(r,:);
end

Answer 3

[Bruno Luong]

[I,J] = ndgrid(1:Nb,1:Nc);
C = accumarray([I(:),J(:)+(I(:)-1)*Nc],A(:));