Parfor大矩阵操作中‘反序列化时内存不足’问题优化求助

问题背景 #

使用parfor执行大型矩阵操作时触发Out of Memory during deserialization错误，同时代码中wim_4d被标记为广播变量，引发不必要的通信开销警告，且无法高效切片wim_4d(:,:,:,ii)。

核心优化方案 #

一、解决内存不足问题 #

1. 放弃4D矩阵转换，保留cell数组直接访问
   原代码中cell2mat将所有图像合并为4D数组，会占用大量连续内存，且导致parfor广播整个大数组。直接保留warpedImages cell数组，在parfor中通过索引获取单张图像：

   % 移除wim_4d = cell2mat(reshape(warpedImages,1,1,1,[]));这一行
   % parfor内直接用warpedImages{ii}替代wim_4d(:,:,:,ii)
   

   此操作让每个worker仅加载所需单张图像，大幅降低内存传输与占用峰值。

2. 重构getIbarNij函数，减少临时大数组
   原函数中创建与图像同尺寸的Ibarij、Ibarji数组属于冗余操作，直接提取重叠区域计算均值即可：

   function [Ibarij,Ibarji,Nij] = getIbarNij(Imij, Imji)
       maski = imbinarize(rgb2gray(Imij));  
       maskj = imbinarize(rgb2gray(Imji));
       
       Nij_im = maski & maskj;
       Nij_im = imfill(Nij_im, 'holes');  
       Nij = sum(Nij_im(:));
       
       if Nij == 0
           Ibarij = zeros(1,3);
           Ibarji = zeros(1,3);
           return;
       end
       
       Nijidx = repmat(Nij_im, 1,1,3);
       % 直接计算重叠区域均值，无需创建全尺寸临时数组
       Ibarij = mean(Imij(Nijidx), 1);
       Ibarji = mean(Imji(Nijidx), 1);
       
       Ibarij(isnan(Ibarij)) = 0;
       Ibarji(isnan(Ibarji)) = 0;
   end
   

   同时提前判断Nij=0的情况，避免无效计算与NaN处理。

3. 用parallel.pool.Constant封装静态数据
   将warpedImages封装为并行池常量，确保每个worker仅加载一次数据，避免重复广播：

   warpedImages_const = parallel.pool.Constant(warpedImages);
   % parfor内访问方式：warpedImages_const.Value{ii}
   

二、消除广播变量警告与切片问题 #

1. 直接用cell数组索引替代4D矩阵切片
   替换原代码中所有wim_4d(:,:,:,x)为warpedImages{x}，彻底消除广播变量警告，同时解决切片效率问题：

   parfor i = 1:length(IuppeIdx)
       [ii,jj] = ind2sub(matSize, IuppeIdx(i));                 
       if ii == jj
           diag_val_1 = 0;
           diag_val_2 = 0;
           Z = setdiff(1:n, ii); % 用setdiff替代手动删除元素，更高效
           for d = Z
               [Ibarij, Ibarji, Nij] = getIbarNij(warpedImages{ii}, warpedImages{d});
               diag_val_1 = diag_val_1 + 2*Nij * Ibarij.^2;
               diag_val_2 = diag_val_2 + Nij;
           end
           diag_val = diag_val_1 + (sigmaN^2/sigmag^2) * diag_val_2;
           B_val = (sigmaN^2/sigmag^2) * diag_val_2;
           Amat_temp{i} = diag_val;
           Bvec(i) = B_val;
       end       
       if ii ~= jj
           [Ibarij,Ibarji,Nij] = getIbarNij(warpedImages{ii}, warpedImages{jj});
           Amat_temp{i} = -2*Nij * (Ibarij .* Ibarji);
       end
   end
   

2. 提前生成上三角索引配对，减少计算开销
   原代码中每次迭代调用ind2sub效率较低，可提前生成所有上三角的ii和jj配对：

   % 替换IuppeIdx的生成逻辑
   [ii_list, jj_list] = find(triu(ones(n,n)));
   % parfor遍历改为：
   parfor i = 1:length(ii_list)
       ii = ii_list(i);
       jj = jj_list(i);
       % 后续逻辑不变
   end
   

三、额外内存优化建议 #

• 及时清理临时变量：在循环内用clear释放不再使用的变量（如Z、maski），减少内存占用。
• 改用单精度计算：若精度允许，将double类型替换为single，可减少一半内存占用。
• 调整worker数量：减少并行worker数量，降低整体内存竞争，避免内存耗尽。

内容的提问来源于stack exchange，提问作者PManjunatha