out-of-memory because of large array: tall array as a workaround?

Given a symmetric, positive definite 3x3 matrix (6 independent components), I want to figure out the spread of the first and second principal invariant.

To this end, I produced the following code and used parfor to speed up the calculations.

%input vectors
nPts = 10;
ub = 3.0;
a = linspace(0.1, ub, nPts); %diagonal entries are > 0
b = linspace(0.1, ub, nPts);
c = linspace(0.1, ub, nPts);
d = linspace(-0.3, ub, nPts); %off-diagonal entries not necessarily > 0
e = linspace(-0.3, ub, nPts);
f = linspace(-0.3, ub, nPts);
%build all combinations
comb = combinations(a,b,c,d,e,f).Variables;
%results array
res = zeros(length(a)*length(b)*length(c)*length(d)*length(e)*length(f), 2);

parfor idx=1:size(comb,1)
    C = [comb(idx,1) comb(idx,4) comb(idx,5); 
        comb(idx,4) comb(idx,2) comb(idx,6); 
        comb(idx,5) comb(idx,6) comb(idx,3)]; %this is the symmetric matrix
    J = sqrt(det(C));
    if J<0.01 
        res(idx, :) = [NaN, NaN];
    else
        Cbar = J^(-2/3).*C; %scaled version of C with det(Cbar)=1
        res(idx, :) =  [trace(Cbar), trace(inv(Cbar))];
    end
end

Starting parallel pool (parpool) using the 'Processes' profile ...
Parallel pool using the 'Processes' profile is shutting down.

%remove rows with NaN (detF<0)
res(isnan(res(:,1)), :) = [];
%remove rows with unphysical behavior
res(res(:,1)<0, :) = [];
res(res(:,2)<0, :) = [];

%create scatter plot
scatter(res(:,1), res(:,2));

This code works, but I have to consider a wider range of the coefficients (a,b,c,d,e,f). For instance, by setting nPts=30 and ub=3.0. If I do this, I ran out-of-memory on my local machine.

Based on the decision table here, using tall arrays may be a workaround in my case. I have all matlab toolboxes installed.

Intuitively, I think my problem can be treated with a tall array. In fact, I do not need to store the entire comb array since my parfor loop operates on each row individually. But I do not know how to translate this into code, if possible at all in my case.

Any help or alternatives are greatly appreciated!

Answer 1

[ProblemSolver]

Hello @SA-W --

You could use tall arrays that you suggested, and then you need to divide the data and then recollected them, and plot it. Here is an example: using the 30 pts that you requested. This also took around 2239.3985 seconds as per my system configurations.

%%
tic

% Define input parameters
nPts = 30;
ub = 3.0;
a = linspace(0.1, ub, nPts); % diagonal entries are > 0
b = linspace(0.1, ub, nPts);
c = linspace(0.1, ub, nPts);
d = linspace(-0.3, ub, nPts); % off-diagonal entries not necessarily > 0
e = linspace(-0.3, ub, nPts);
f = linspace(-0.3, ub, nPts);

% Create a parallel pool
parpool();

% Initialize variables
chunkSize = 10000; % Adjust the chunk size as per your memory capacity
totalCombinations = numel(a) * numel(b) * numel(c) * numel(d) * numel(e) * numel(f);
numChunks = ceil(totalCombinations / chunkSize);
resCell = cell(numChunks, 1);
chunkIdx = 1;

% Process combinations in smaller chunks
for aIdx = 1:numel(a)
    for bIdx = 1:numel(b)
        for cIdx = 1:numel(c)
            for dIdx = 1:numel(d)
                for eIdx = 1:numel(e)
                    for fIdx = 1:numel(f)
                        % Construct symmetric matrix
                        C = [a(aIdx) d(dIdx) e(eIdx);
                             d(dIdx) b(bIdx) f(fIdx);
                             e(eIdx) f(fIdx) c(cIdx)];
                        J = sqrt(det(C));
                        if J >= 0.01
                            Cbar = J^(-2/3) * C;
                            invValues = [trace(Cbar), trace(inv(Cbar))];
                            resCell{chunkIdx} = [resCell{chunkIdx}; invValues];
                        end
                        
                        % Move to the next chunk if necessary
                        if size(resCell{chunkIdx}, 1) >= chunkSize
                            chunkIdx = chunkIdx + 1;
                        end
                    end
                end
            end
        end
    end
end

% Concatenate the results from all chunks
res = cell2mat(resCell);

% Remove rows with NaN (detF < 0)
res(isnan(res(:,1)), :) = [];

% Remove rows with unphysical behavior
res(res(:,1) < 0, :) = [];
res(res(:,2) < 0, :) = [];

% Create scatter plot
scatter(res(:,1), res(:,2));

% Delete the parallel pool
delete(gcp);

toc

[SA-W]

Thank you for your answer, however I have some queries:

% Concatenate the results from all chunks
res = cell2mat(resCell);

You create here a variable res which stores the results of all combinations. Does this not result into the out-of-memory error that I wanted to avoid for a large value of nPts?

What happens if you set nPts=60 ? If I am not mistaken, res would require roughly 60^6*2*8*1e-9 = 746 GB memory. Currently, we only require around 11 GB memory, which your operating system can give you.

[ProblemSolver]

@SA-W -- You are correct. What you can do for that is process the results in chunck and write those chuncks into the tempDir. In this way you can limit the memory usage to the current chunk being processed. Something like this:

% Create a temporary directory for intermediate results
tempDir = fullfile(tempdir, 'intermediate_results');
if exist(tempDir, 'dir') ~= 7
    mkdir(tempDir);
end

Then within the loop you can save the result. You will need to create a processing chunk code where, it will handle each chunk of the saved result. Once that chunk size is reached, the function or the section of the code will save it and concatenates it and stores it in a variable res.

You have to preallocate main result array something like this:

res = [];
save('result.mat', 'res');

In this way you can keep increasing the points, and still you will be able to perform the operation that you need.

[SA-W]

Once that chunk size is reached, the function or the section of the code will save it and concatenates it and stores it in a variable res.

But if I concatenate the individual chunks, I end up again with the large array, right? Would you please extend your answer by the save operations to make the overall strategy more clear for me?

[ProblemSolver]

But if I concatenate the individual chunks, I end up again with the large array, right?

Yes, if you concatinate the individual chunks it will still create a large array, but then it will depend upon the your HARDDRIVE space and not the RAM memory (the error you get "out-of-memory").

Would you please extend your answer by the save operations to make the overall strategy more clear for me?

As you pointed out that the concatenation of all results into a single array 'res' can indeed result in an out-of-memory error for very large values of 'nPts'. The memory required to store this said results array would exceed the available memory capacity. To avoid this issue, I provided you an approach that processes the results in chunks (meaning it takes say 200 points out of the one large array) and write that chunk into the HARDDRIVE space before moving or processing the next chunk (i.e., say another 200 points out of the same large array).

The temperary directory part that I suggested:

What it will do is when the function takes this said 200 points, it stores these intermediate results in a separate ".MAT" file. In layman terms: "instead of storing all the intermediate results in memory, the code saves each chunk of intermediate results to separate files on your computer's storage. This way, the memory usage remains low because only one chunk of results is loaded and processed at a time. Once a chunk is processed and its results are incorporated into the final output, the temporary file is deleted to free up disk space. By saving intermediate results to disk, the code can handle larger amounts of data that would otherwise exceed the available memory. It allows the computation to be performed in smaller manageable chunks, reducing the overall memory requirements and ensuring that the code doesn't run out of memory. Overall, saving intermediate results to disk is a strategy to handle large amounts of data by temporarily storing and retrieving partial results from storage, thus enabling the code to handle larger computations without running into memory limitations."

Now the second of this approach is loading these chucks and processing these chucks:

the intermediate results are saved as separate MAT files for each chunk (as mentioned in the above point). Once the parallel processing of all combinations is complete and the parallel pool is deleted, the code moves on to the step of loading and processing each chunk individually. The main purpose of this step is to retrieve the intermediate results from disk and incorporate them into the final result. By loading and processing one chunk at a time, the code ensures that the memory usage is limited to the current chunk being processed, preventing the accumulation of large arrays in memory. The approach that I provided iterates over each chunk, which is represented by a unique MAT file. For each iteration, you will construct the full path of the MAT file and checks if it exists. If the MAT file exists, it will be loaded using the load function, and the relevant data will be extracted and appended to the res array (that is your final array for scatter plot). After processing each chunk, the corresponding MAT file should be deleted to free up disk space.

So in other words: loading and processing each chunk separately is a memory-efficient approach that will allow your code to handle large datasets by loading and processing a portion of the data at a time, reducing memory usage and preventing memory-related issues.

I hope this explanation helps!

[SA-W]

Thanks for your comphrehensive explanation. Let me ask a follow-up regarding

If the MAT file exists, it will be loaded using the load function, and the relevant data will be extracted and appended to the res array

The overall strategy we want to pursue is clear to me, but when we loop over all .MAT files and append their contents to the final "res" array, we have again the large "res" array at the end of the loop, do not we?

[ProblemSolver]

Yes. Thats how it will plot the graph.

[SA-W]

But then it can not work if we assume that the "res" array does not fit into memory? We basically shifted the problem after the processing step right?

[ProblemSolver]

I guess you don't care about the intermediate results for additional analyses, or computation or anything.

You can just plot simulataneously as your process your data. Then your code is already loop on each row, you can just keep appending the scatter plot as your process. I don't recommend this as this will hog your system. But I guess that is what you want.

% Update the scatter plot with current results
hold on;
scatter(invValues(1), invValues(2), 'b.');
hold off;
drawnow;

I hope this helps!