Effacer les filtres
Effacer les filtres

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

4 vues (au cours des 30 derniers jours)
SA-W
SA-W le 12 Juil 2023
Commenté : ProblemSolver le 14 Juil 2023
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!

Réponse acceptée

ProblemSolver
ProblemSolver le 12 Juil 2023
Modifié(e) : ProblemSolver le 12 Juil 2023
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
  8 commentaires
SA-W
SA-W le 13 Juil 2023
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
ProblemSolver le 14 Juil 2023
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!

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Matrix Indexing dans Help Center et File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by