Out of memory - an alternative to this algorithm?

I am working on Markov Chain problem involving a large, sparse transition matrix. Unfortunately, in generating the transition matrix as I do, my system quickly runs out of memory. Is my algorithm (seen below) causing this problem unnecessarily early, or is it just the nature of the problem that prevents computation on my system?

n = 5;
L = 0.1;
r = [];
c = [];
s = [];
for i = 1:combinations % combinations the number of combinations in E_perms, a matrix holding (1/L + 1)^n possible states. This can obviously be extremenly large.
     state_0 = u(i,:);
    [states_1,transition_probs] = transitionProbability(state_0,k,L,E_perms); % A function that generates a vector of possible states after state_0, and the associated probabilites of transitioning to them. Note that only a small subset of all states can be reached from state_0 so the final transition matrix is very sparse. 
    for ii = 1:size(states_1(:,1))
        [~,indx] = ismembertol(states_1(ii,:),u,0.00001,'ByRows', true);
        r = [r i];
        c = [c indx];
        s = [s transition_probs(ii,1)];
    end
end
transition_matrix = accumarray([r',c'],s',[],[],[],1);

[KSSV]

combinations is not given......a function is also not given..

Answer 1

[Walter Roberson]

Your indx values could be up to numel(u) and your ii values can be up to combinations. That should give you enough information to be able to spalloc() a result matrix before doing any looping.

Now, as you go through, do not build a overall r, c, s matrices: build them per i. Then after the ii loop, accumarray what you got from that loop, and specify the known maximum sizes as your third parameter for accumarray . Continue to specify sparse for that matrix. After you have built that matrix with the same sparse bounds as the transition_matrix , add it to the transition_matrix .

Your working memory would then be only enough for one generation, in linear form, plus the sparse memory for representing that linear in sparse form, plus the memory for the sparse final matrix. I guess potentially you could need up to twice the memory of the sparse final matrix as it is probably not going to do the sparse addition "in place".

Your current system requires three arrays full of linear information, plus enough memory to store the sparse final array, so by my calculation building iteratively should require less memory.

[Ulrik William Nash]

Thank you Walter, I see what you mean. I am currently trying to implement this part of you advice:

"...do not build a overall r, c, s matrices: build them per i. Then after the ii loop, accumarray what you got from that loop".

I see the effect of this, but I am unsure how to do a forth-running accumarray. Could you please show how that is done? (I understand the intuition: In my code, I create long vectors and then accumarray at the very end. What you advise me to do is compress at each iteration.)

[Walter Roberson]

You might have enough information to give a better formula for the total number of non-zero entries that you will end up with; put that in nzguess

maxrow = combinations;
maxcol = numel(u);
nzguess = floor(maxrow * maxcol / 1000);

transition_matrix = spalloc(maxrow, maxcol, nzguess);

    n = 5;
    L = 0.1;

for i = 1:combinations % combinations the number of combinations in E_perms, a matrix holding (1/L + 1)^n possible states. This can obviously be extremenly large.
    r = [];
    c = [];
    s = [];
    state_0 = u(i,:);
    [states_1, transition_probs] = transitionProbability(state_0,k,L,E_perms); % A function that generates a vector of possible states after state_0, and the associated probabilites of transitioning to them. Note that only a small subset of all states can be reached from state_0 so the final transition matrix is very sparse. 
    for ii = 1:size(states_1(:,1))
        [~,indx] = ismembertol(states_1(ii,:),u,0.00001,'ByRows', true);
        r = [r i];
        c = [c indx];
        s = [s transition_probs(ii,1)];
    end
    transition_submatrix = accumarray([r',c'], s', [maxrow, maxcol], [], [], 1);
    transition_matrix = transition_matrix + transition_submatrix;
end

[Ulrik William Nash]

By the way, how important is an accurate allocation of nnz? And is overshooting worse than undershooting?

[Walter Roberson]

You posted a Question about that; James answered (correctly) that undershooting can be much worse than overshooting.

[Ulrik William Nash]

Thank you Walter, yes, I did.

Unfortunately, after implementing your solution of "chopping and compressing" more frequently than I had done, the memory consumption of MATLAB is surprisingly out of tune with the size of the matrices being generated. This I cannot understand (see attached image, which shows the memory usage climbing from the start of the script till it drops at the finish).

Do you have any idea what might be going on? And any idea about solutions? I have thought about crunching segments of states, saving the output, moving to the next segment of states, and so on, with clearvars in between. Is that wise?

[Walter Roberson]

Be cautious: clearing variables in a loop is expensive. I suspect is less expensive to assign empty array to them, but that would need to be timed.

Answer 2

[Jan]

This note seems to be important:

combinations the number of combinations in E_perms, a matrix holding
(1/L + 1)^n possible states. This can obviously be extremenly large.

For the given data, n = 5, L = 0.1, you get 161'051. Is this "extremely large" already? How large is E_perms and u? What is k? Is transitionProbability a user defined function?