2

In MATLAB have a large matrix with transition probabilities transition_probs, and an adjacency matrix adj_mat. I want to compute the cumulative sum of the transition matrix along the columns and then element wise multiply it against the adjacency matrix which acts as a mask in this way:

cumsumTransitionMat = cumsum(transition_probs,2) .* adj_mat;

I get a MEMORY error because with the cumsum all the entries of the matrix are then non-zero.

I would like to avoid this problem by only having the cumulative sum entries where there are non zero entries in the first place. How can this be done without the use of a for loop?

Vass
  • 2,682
  • 13
  • 41
  • 60

2 Answers2

2

You can apply cumsum on the non-zero elements only. Here is some code:

A = sparse(round(rand(100,1)));        %some sparse data
A_cum = A;                             %instantiate A_cum by copy A

idx_A = find(A);                       %find non-zeros
A_cum(idx_A) = cumsum(A(idx_A)); %cumsum on non-zeros elements only

You can check the output with

 B = cumsum(A);


   A_cum   B
     1     1
     0     1
     0     1
     2     2
     3     3
     4     4
     5     5
     0     5
     0     5
     6     6

and isequal(A_cum(find(A_cum)), B(find(A_cum))) gives 1.

marsei
  • 7,691
  • 3
  • 32
  • 41
  • I am not sure if this is working. The results are different using A_cum – Vass Sep 13 '13 at 08:11
  • 1
    Please see the edit - results are identical to me: non-zeros A elements and their corresponding B are identical. Unless I didn't correctly understand either your comment or question. Could you elaborate? – marsei Sep 13 '13 at 08:23
2

when CUMSUM is applied on rows, for each row it will go and fill with values starting with the first nonzero column it finds up until the last column, thats what it does by definition.

The worst case in terms of storage is when the sparse matrix contains values at the first column, the best case is when all nonzero values occur at the last column. Example:

% worst case
>> M = sparse([ones(5,1) zeros(5,4)]);
>> MM = cumsum(M,2);       % completely dense matrix
>> nnz(MM)
ans =
    25

% best case
>> MM = cumsum(fliplr(M),2);

If the resulting matrix does not fit in memory, I dont see what else you can do, except maybe use a for-loop over the rows, and process the matrix is smaller batches...

Note that you cannot apply the masking operation before computing the cumulative sum, since this will alter the results. So you cant say cumsum(transition_probs .* adj_mat, 2).

Amro
  • 123,847
  • 25
  • 243
  • 454
  • so it's not possible? do you know what Im trying to do? I want to make a Markov transition matrix. So in a row (state) there are only a certain set of non-zero entries ie, possible state transitions denoted by a adjacency matrix of 1/0. I have a set of probabilities for those state transitions. I throw a random number for a candidate next state and want to see which column number values it falls between in the cumulative sum. How else can I do this? – Vass Sep 13 '13 at 21:51
  • well it seems like I perform the operation row by row calculating it each time there is a transition – Vass Sep 17 '13 at 13:26
  • yes I think I understood what the code is supposed to do (basically you want to pick a next state according to a vector of weights corresponding to transition probabilities of the current state, but restricted to the adjacent states only). You are sort of implementing [`randsample`](http://stackoverflow.com/q/2977497/97160). The problem is that when the result is computed over the entire matrix at once, it is too dense to fit in memory. That's why I suggested looping over the rows one by one instead of a vectorized call.. Sorry I couldn't be of more help :) – Amro Sep 17 '13 at 13:56
  • 1
    on second thought, I think there might be room to improve the code. But first can you explain in more details what you are doing? I dont get why you are computing the cumulative sum for the entire matrix in the first place.. If you know you are in state `i` you only need to look at its corresponding row `transition_probs(i,:)`, is my understanding correct? – Amro Sep 17 '13 at 14:18
  • 1
    yes that is correct, I now have the code `cumsumTransitionRow = cumsum(transitionMat(row_temp,:)) . graph_com(row_temp,:); nodesLess = find(rand <= cumsumTransitionRow(:));` so I compute the row of interest for that state `i`, the appropriate row number of the transition matrix, mask it and use the commulative sum of the transition probabilities to sample the next state – Vass Sep 17 '13 at 15:41