Inserting zero rows and columns into matrix efficiently

I'm currently working on a project where I have a large set of matrices of varying size. Each matrix, A, has a vector, pointer, which corresponds to location information. I'm trying to get all the matrices into the same "space"; to achieve this, I am looking to create a "master" matrix, M, which has all the elements of A expanded to the larger matrix, which has location data for all matrices in the set. I have already created a master pointer, pointer_master, which contains the locations found in all matrices.

I'm trying to create two matrices using the original matrix, A, pointer, and pointer_master. The first matrix is the "occurrence" matrix, O, which shows you what elements from A appear in the large matrix M. I have already done this by using the function setdiff to find the missing elements

Code:

 % This code compares pointer_master to pointer and creates
 % a logical matrix O which sets all missing indices to 0.
 pointer_master = 1:25;
 pointer = [1 4 8 10 11 12 13 14 15 17 21 23 25];
 [~, pointer_missingIndex] = setdiff(pointer_master,pointer);
 O = true(length(pointer_master));
 O(:,pointer_missingIndex) = false;
 O(pointer_missingIndex,:) = false;

Here we have a pointer that is missing 12 elements from the master pointer. The code above creates a new large matrix O, with 0 in place of these missing elements. When this process is complete, I intend to sum these matrices to get the frequency of each element (I also want to average across elements using this matrix).

I'm having trouble figuring out how to create the master matrix. Unlike the occurence matrix, I need the elements of A to appear in M with all its elements. I can do this inefficiently, but I'm sure there is a better way to do this than what I have right now.

Code:

 M = A;
 for ii = 1:length(pointer_master)
   X = find(pointer==pointer_master(ii), 1);
   if(isempty(X))
      column_insert = zeros(size(M,2),1); %#ok<NASGU>
      row_insert = zeros(1,size(M,2)+1); %#ok<NASGU>
      if(pointer_master(ii)<min(pointer))
         column_insert = zeros(size(M,2),1);
         row_insert = zeros(1,size(M,2)+1);
         M = [logical(column_insert) M]; %#ok<AGROW>
         M = [logical(row_insert); M]; %#ok<AGROW>
       elseif(pointer_master(ii)>max(pointer))
         column_insert = zeros(size(M,2),1);
         row_insert = zeros(1,size(M,2)+1);
         M = [logical(column_insert) M]; %#ok<AGROW>
         M = [logical(row_insert); M]; %#ok<AGROW>
       else
         column_insert = zeros(size(M,2),1);
         row_insert = zeros(1,size(M,2)+1);
         M = [M(:,1:N) logical(column_insert) M(:,N+1:end)];
         M = [M(1:N,:); logical(row_insert); M(N+1:end,:)];
       end
     else
       N = ii;
     end
end

I know this method is very inefficient, especially for matrices with 20,000+ elements, so I'm trying to figure out how to do this faster. I can use ismember to output a vector that tells you where pointer falls in pointer_master. but I don't know how to get the elements from A into M.

Code:

 [~,X] = ismember(pointer_master,pointer);

Does anyone have any experience with this? Thanks.

For example:

 A = [1 2 3; 4 5 6; 7 8 9]; 
 pointer = [1 3 4]; 
 pointer_master = 1:5;

My occurrence matrix would be:

 O = [1 0 1 1 0; 0 0 0 0 0; 1 0 1 1 0; 1 0 1 1 0; 0 0 0 0 0];

My master matrix M would be:

 M = [1 0 2 3 0; 0 0 0 0 0; 4 0 5 6 0; 7 0 8 9 0; 0 0 0 0 0];

[Sean de Wolski]

Could you perhaps boil that down to a (small example)?

It's not clear why you could not preallocate an M matrix and then just set it as necessary.

[Sean de Wolski]

Also, what is A?

[Effrom]

I added an example to my original question, *A*, is an input matrix that has non-zero elements.

As for master matrix generation code, I will likely not use it all, I just wrote code that would work. I would have pre-allocated the matrix, but I'm not sure how to insert the zero rows and columns beforehand.

Answer 1

[Andrei Bobrov]

A = [1 2 3; 4 5 6; 7 8 9]; 
pointer = [1 3 4]; 
pointer_master = 1:5;

p = ismember(pointer_master,pointer);
M = bsxfun(@and,p,p')+0;
M(M~=0)=A;

Answer 2

[Richard Brown]

Perhaps I'm misunderstanding, but your example can be generated from the following

M = zeros(numel(pointer_master));
M(pointer,pointer) = A;

If the number of elements in pointer is much fewer than the number of elements of pointer_master, you may well want to use a sparse matrix:

n = numel(pointer_master);
[I, J] = ndgrid(pointer);
M = sparse(I, J, A, n, n);

[Andrei Bobrov]

+1

Answer 3

[Jan]

I do not understand what you are looking for. But I have a general advice when I see "%#ok<AGROW>". Don't do this. Either pre-allocate the required size and fill in the data afterwards. Or pre-allocate to much memory and cut off the unneeded part at the end.

This wastes time also:

column_insert = zeros(size(M,2),1);
row_insert = zeros(1,size(M,2)+1);
M = [M(:,1:N) logical(column_insert) M(:,N+1:end)];

At first "false()" is faster than "logical(zeros())". But I assume this would be even faster:

M(:, N+2:end+1) = M(:,N+1:end);
M(:, N+1) = 0;

I have the impression that you mix DOUBLEs and LOGICALs in the matrix M. If so, don't do it.

Answer 4

[Effrom]

Thank you all for your help.

Richard, I initially tried M(pointer,pointer) = A, but I ran into memory issues (I'm working with large 20,000+ element matrices). Thank you Andrei, the bsxfun (I'm always learning about new functions) seems to be the best solution. I also came up with my own solution, that while more efficient than what I started with, is not as efficient as your solution.

My Solution: Since the matrices I deal with contain zero elements, I created a new matrix A:

 A =

     0     2     0   | pointer 1
     4     0     6   |         3
     0     8     0   |         4

The matrix A is paired with a pointer that contains location information. The pointer (pointer = [1 3 4]) corresponds to locations 1, 3 and 4. The master matrix contains the locations 1, 2, 3, 4 and 5 (pointer_master = 1:5). So the goal is to map the original elements of A to new elements in the master matrix M:

 M =

     0     0     2     0     0   | pointer_master 1
     0     0     0     0     0   |                2
     4     0     0     6     0   |                3
     0     0     8     0     0   |                4
     0     0     0     0     0   |                5

To get from A to M, I generated the occurrence matrix O as I described earlier:

 [~, pointer_missingIndex] = setdiff(pointer_master,pointer);
 O = true(length(pointer_master));
 O(:,pointer_missingIndex) = false;
 O(pointer_missingIndex,:) = false;

 O =

     1     0     1     1     0   | pointer_master 1
     0     0     0     0     0   |                2
     1     0     1     1     0   |                3
     1     0     1     1     0   |                4
     0     0     0     0     0   |                5

I now realize the bsxfun function achieves the same outcome. To create the master matrix, I decided to vectorize my matrices:

 A_v = reshape(A,numel(A),1); % Vector of input matrix
 O_v = reshape(O,numel(O),1); % Vector of occurrence matrix
 M_v = zeros(numel(O),1);     % Vector of master matrix (empty)

I then found the non-zero elements in O, which gives you the indices of A as they would appear in M. I wanted to work with the non-zero elements only, so I added these last few steps:

 O_vrow1 = find(O_v==true); % Non-zero elements in O
 A_vrow = 1:numel(A);       % All elements of A
 A_vrow0 = find(A_v==0);    % Zero elements in A
 A_vrow1 = find(A_v>0);     % Non-zero elements in A

 M_vrow1 = O_vrow1;
 M_vrow1(A_vrow0) = [];      % Remove zero elements from M_vrow1
 M_v(M_vrow1) = A(A_vrow1);  % Insert elements of A into M_v
 n = length(pointer_master);
 M = reshape(M_v,n,n);       % Master matrix

The idea behind these last steps was to generate a vector containing the non-zero elements of A in terms of M. I got a solution that is more efficient, but I see even better solutions here. Thanks to all again.

[Richard Brown]

Glad to hear you got it sorted. I'm puzzled as to why the bsxfun version works differently for you - on my system both methods use exactly the same amount of memory (the amount of memory required to store both A and M as doubles).

[Effrom]

I'm not sure why it didn't work because it worked on the smaller matrices, but didn't seem to scale up.