How to update variable within a matfile inside a parfor loop?

I have a very expensive loop that I'm trying to parallelize, and part of this loop involves updating an entry in a 4D array inside a matfile (I must save results to disk and access them through a matfile pointer due to RAM limitations). However, I get an error that says that the matfile pointer variable cannot be classified. As an illustrative example of what I'm trying to do, consider the code below:

testOut = []; % create variable to try and update
save('TestFile.mat', 'testOut', '-v7.3'); % save variable into accessible matfile
FileOut = matfile('TestFile.mat','Writable',true); % Set up pointer to matfile
nrows = 200; ncols = 200; nplanes = 32; nvolumes = 10; % Set up dimensions of 4D array
FileOut.testOut = single(zeros(nrows,ncols,nplanes,nvolumes)); % Set initial size of variable in matfile
parfor i = 1:nvolumes
   FileOut.testOut(:,:,:,i) = i; % artificial example, point is that I want to update variable using fourth dimension index
end

The error this code would report is:

 Error: The variable FileOut in a parfor cannot be classified.
 See Parallel for Loops in MATLAB, "Overview".

Basically, the loop is performing an independent calculation on a 3D volume image each time, and each resulting volume image is saved in a 4D array using the fourth dimension to mark volume image number. However, I often have thousands of these images, and so the 4D array must be saved to disk and accessed from a matfile to avoid overloading the RAM. I'd like to adapt my code to use parfor, but I can't figure out how to get parfor to play nicely with the matfile pointer. Can anyone help me out here, please? I understand that sliced variables must be used with parfor, and that variables of the form I'm using here are not allowed, but I can't figure out a solution...

Answer 1

[Edric Ellis]

It might be possible to overcome the "slicing" problems you're seeing here - but you're still left with the fundamental underlying problem that you're trying to get multiple worker processes to write to the same file concurrently. That is never going to work well, as the writes will conflict and almost inevitably corrupt your file.

What I'd suggest is having each worker write to a temporary .mat file during the parfor loop, and then run a post-processing stage to collect the results. (I'm presuming here that computing stuff to go into the file takes a long time, but accessing the data in the file is relatively inexpensive). I'm going to use the parallel.pool.Constant from R2015b, but the same can be achieved using the Worker Object Wrapper.

This example is a bit involved, hopefully it's clear what's going on. You'll need to adapt things a little to get them to work with your multi-dimensional data.

%%step 1: create a mat-file per worker using SPMD
spmd
    myFname = tempname(); % each worker gets a unique filename
    myMatfile = matfile(myFname, 'Writable', true);
end

%%step 2: create a parallel.pool.Constant from the 'Composite'
% This allows the worker-local variable to used inside PARFOR
myMatfileConstant = parallel.pool.Constant(myMatfile);

%%Step 3: run PARFOR
parfor idx = 1:100
    resultToSave = idx * 100;
    matfileObj = myMatfileConstant.Value;
    % Append into 'testOut', storing the index
    matfileObj.testOut(1, idx) = resultToSave;
    matfileObj.gotResult(1, idx) = true;
end

%%Step 4: accumulate the results on the client
% Here we retrieve the filenames from 'myFname' Composite,
% and use them to accumulate the overall result
outmatfile = matfile('out.mat', 'Writable', true);
for idx = 1:numel(myFname)
    workerFname = myFname{idx};
    workerMatfile = matfile(workerFname);
    workerOutSz = size(workerMatfile, 'testOut');
    for jdx = 1:workerOutSz(2)
        if workerMatfile.gotResult(1, jdx)
            outmatfile.out(1, jdx) = workerMatfile.testOut(1, jdx);
        end
    end
end

[Nathan Murtha]

Thanks Edric, great answer. I thought something like this might end up being my only option, thanks very much for your detailed response. I implemented this, adapting it for my 4D array, and it executed just fine. However, it seems like the results differ from those produced from my serial code! To put it simply, the code is supposed to reconstruct a 3D volume image, storing them all in a 4D array. When I run the serial and parallel versions of the code, the results differ slightly (which I can verify by taking the difference of each individual image). Any idea why this might be? Would there be some kind of error introduced through the handling of numbers by each CPU? I wish that I could ask a better question here, and I can certainly provide more details if that would help.

[Edric Ellis]

The workers executing your parfor loop body differ only slightly from the MATLAB client. Workers run in single-threaded mode - this can introduce slight numerical variations. Workers use different random number generators, this can introduce more noticeable discrepancies.

[Timur Mokaev]

Dear Edric. Could you clarify how to modify this example when my 'resultToSave' is not a single value but a vector of size [1, N] or a matrix of size [M, N]. Also, I can't understand when you determine 'testOut' and 'gotResult' fields for a mat objects at each worker. Within the spmd we just create mat files without any data.

[Edric Ellis]

Here's an updated version of the answer. I'm storing the data in a cell array, and I've modified things to pre-allocate the variables in the matfile objects. I've also changed things slightly to accumulate the final result in memory. But you could equally assign directly into another matfile for each page.

%%step 1: create a mat-file per worker using SPMD
N = 100;
spmd
    myFname = tempname(); % each worker gets a unique filename
    myMatfile = matfile(myFname, 'Writable', true);
    % Seed the variables in the matfile object
    myMatfile.testOut = cell(1, N);
    myMatfile.gotResult = false(1, N);
end

%%step 2: create a parallel.pool.Constant from the 'Composite'
% This allows the worker-local variable to used inside PARFOR
myMatfileConstant = parallel.pool.Constant(myMatfile);

%%Step 3: run PARFOR
m = 3;
n = 2;
parfor idx = 1:N
    resultToSave = rand(m, n);
    matfileObj = myMatfileConstant.Value;
    % Assign into 'testOut'
    matfileObj.testOut(1, idx) = {resultToSave}
    matfileObj.gotResult(1, idx) = true;
end

%%Step 4: accumulate the results on the client
% Here we retrieve the filenames from 'myFname' Composite,
% and use them to accumulate the overall result in memory

% We're going to concatenate results in the 3rd dimension
result = nan(m, n, N);
for idx = 1:numel(myFname)
    workerFname = myFname{idx};
    workerMatfile = matfile(workerFname);
    for jdx = 1:N
        if workerMatfile.gotResult(1, jdx)
            pageCell = workerMatfile.testOut(1, jdx);
            result(:, :, jdx) = pageCell{1};
        end
    end
end

[Timur Mokaev]

Thanks a lot!

[Aditya Nanda]

Hi Edric! this is a great answer. I have a very similar issue that I am struggling with-> the only difference being my "resulttoSave" is multidimenisonal

For instance, here is the code (this is Step 3 in the above example)

parfor idx = 1:100
    resultToSave = someFunction; 
    matfileObj = myMatfileConstant.Value;
    % Append into 'testOut', storing the index
    matfileObj.testOut(1:n,1:t, idx) = resultToSave;
    matfileObj.gotResult(1, idx) = true;
end
% % 

The variable resultToSave is 3 dimensional (double ). When I run it like this, I get the error :

Variable resultToSave has 2 dimensions thus indexing in third dimension is not possbile. Please help. I am not experienced with using spmd. Thanks

[Edric Ellis]

Hi @Aditya Nanda, please could you post a new question with some slightly more detailed reproduction steps that reproduce the problem in a standalone way. (Feel free to pop a link in a comment here so I'll get notified).

Answer 2

[Jason Climer]

It's worth noting that accumulating the results into a matfile instead of the local memory is prohibitively slow.

[Timur Mokaev]

Sure, but I guess, here we consider the case when accumulated results may not fit in local memory. Also, if for some reason the whole numerical procedure crashes, one will lose all the results stored in local memory.