加速巨大稀疏矩阵的条件填充

Question

我想知道是否有一种方法可以加速（也许通过矢量化？）巨大稀疏矩阵（例如〜1e10 x 1e10）的条件填充。这是示例代码，其中我有一个嵌套循环，并且仅在满足特定条件时才填充稀疏矩阵：

% We are given the following cell arrays of the same size:
% all_arrays_1
% all_arrays_2
% all_mapping_arrays

N = 1e10; 

% The number of nnz (non-zeros) is unknown until the loop finishes
huge_sparse_matrix = sparse([],[],[],N,N);

n_iterations = numel(all_arrays_1);

for iteration=1:n_iterations

  array_1 = all_arrays_1{iteration};
  array_2 = all_arrays_2{iteration};

  mapping_array = all_mapping_arrays{iteration};

  n_elements_in_array_1 = numel(array_1);
  n_elements_in_array_2 = numel(array_2);

  for element_1 = 1:n_elements_in_array_1

    element_2     = mapping_array(element_1);

    % Sanity check:
    if element_2 <= n_elements_in_array_2
       item_1 = array_1(element_1);
       item_2 = array_2(element_2);      

       huge_sparse_matrix(item_1,item_2) = 1;
    end     
  end   
end

我正在努力对嵌套循环进行矢量化。据我了解，当要填充的条目数量很大（~100M）时，逐个元素填充稀疏矩阵的速度非常慢。我需要使用稀疏矩阵，因为它的尺寸在 10,000M x 10,000M 范围内。然而，这种在 MATLAB 中填充稀疏矩阵的方法非常慢。

编辑：

我更新了变量的名称，以更好地反映它们的性质。没有函数调用。

附录：

此代码为一个巨大的图构建矩阵邻接。变量all_mapping_arrays以本地表示的方式保存图形节点之间的映射数组（〜邻接关系），这就是为什么我需要array_1和< code>array_2 将邻接映射到全局表示。

Answer 1

我认为这将是稀疏矩阵的增量更新，而不是基于循环的条件，这将减慢速度。

当您通过 A(i,j) = 1 之类的方式向稀疏矩阵添加新条目时，通常需要重新打包整个矩阵数据结构。这是一项昂贵的操作。如果您感兴趣，MATLAB 在内部使用 CCS 数据结构（压缩列存储），这在数据结构部分 此处。注意声明：

这种方案对于一次操作一个元素的矩阵来说效率不高。
时间

，最好（更快）地将矩阵中的非零条目累积为一组三元组，然后对 sparse 进行一次调用。例如（警告 - 大脑编译的代码！）：

% Inputs:
% N 
% prev_array and  next_array
% n_labels_prev and n_labels_next
% mapping

% allocate space for matrix entries as a set of "triplets"
ii = zeros(N,1);
jj = zeros(N,1);
xx = zeros(N,1);
nn = 0;

for next_label_ix = 1:n_labels_next

    prev_label     = mapping(next_label_ix);

    if prev_label <= n_labels_prev
        prev_global_label = prev_array(prev_label);
        next_global_label = next_array(next_label_ix);   

        % reallocate triplets on demand
        if (nn + 1 > length(ii))
            ii = [ii; zeros(N,1)];
            jj = [jj; zeros(N,1)];
            xx = [xx; zeros(N,1)];
        end   

        % append a new triplet and increment counter
        ii(nn + 1) = next_global_label; % row index
        jj(nn + 1) = prev_global_label; % col index
        xx(nn + 1) = 1.0;               % coefficient
        nn = nn + 1;
    end     
end

% we may have over-alloacted our triplets, so trim the arrays
% based on our final counter
ii = ii(1:nn);
jj = jj(1:nn);
xx = xx(1:nn);

% just make a single call to "sparse" to pack the triplet data
% as a sparse matrix object
sp_graph_adj_global = sparse(ii,jj,xx,N,N); 

我一次分配 N 条目的块。假设您对矩阵的结构了解很多，您可能可以在此处使用更好的值。

希望这有帮助。