通过 3D 数组页面的交错复制构建邻接矩阵

Question

背景

我正在尝试建模一个可以在每个时间步更改其配置的系统。配置的多样性是预先知道的，不依赖于时间步长。在某些配置之间允许转换，在其他配置之间禁止转换。 objective 是构建允许转换的邻接连接矩阵，它跨越多个时间步长。

设置

设A为表示允许转换的s*s*k逻辑矩阵，A1...Ak表示A的pages/slices:

A1 = A(:,:,1); A2 = A(:,:,2); ... Ak = A(:,:,k);

第3个维度的意思是一个转换需要多少时间步，例如：如果A(1,3,2)不为零，则表示状态#1可以转换到状态#3，这将需要 2 个时间步。

让B成为我们要构建的邻接矩阵，它代表nt个时间步长。 B的形状应该是示意性的（分块矩阵表示法）：

     _                                   _
    | [0] [A1] [A2] ... [Ak] [0]  ... [0] |
B = | [0] [0]  [A1] [A2] ... [Ak] ... [0] |
    |  ⋮    ⋮     ⋱    ⋱      ⋱       ⋮  |
    |_[0] [0]  …  …  …  …  …  …  …  … [0]_| "[A1] [A2] ... [Ak]"

其中主块对角线由 nt 个 0 块组成，A 的切片逐渐向右“推”到“时间用完”，[= 的切片13=] 结束于 B 的“外部”⇒ 表示不再可能进行转换。由于 B 由 nt*nt s*s 个块组成，因此其大小为 (nt*s)×(nt*s).

Question: Given A and nt, how can we construct B in the most CPU- and memory-efficient way?

备注

由于 B 大部分由零填充，因此 sparse.
CPU 在我的应用程序中效率（运行时）比内存效率更重要。
在真题中，s=250和nt=6000。
欢迎外部 scripts/classes/tools。
我的想法不是最初构建矩阵交错，而是有一个 [A1] 块的主对角线和 [=40=]-ing 和掩码，当其他一切都完成时。

演示 + 简单的实现

s = 3; k = 4; nt = 8;
A = logical(cat(3, triu(ones(s)), eye(s), zeros(s), [0 0 0; 0 0 0; 0 1 0]));
% Unwrap A (reshape into 2D):
Auw = reshape(A, s, []);
% Preallocate a somewhat larger B:
B = false(nt*s, (nt+k)*s);
% Assign Auw into B in a staggered fashion:
for it = 1:nt
  B( (it-1)*s+1:it*s, it*s+1:(it+k)*s ) = Auw;
end
% Truncate the extra elements of B (from the right)
B = B(1:nt*s, 1:nt*s);
spy(B);

导致：

Answer 1

一种解决方案是使用隐式扩展计算所有索引：

% Dev-iL minimal example
s = 3; k = 4; nt = 8;
A = logical(cat(3, triu(ones(s)), eye(s), zeros(s), [0 0 0; 0 0 0; 0 1 0]));
Auw = reshape(A, s, []);

% Compute the indice:
[x,y] = find(Auw);
x = reshape(x+[0:s:s*(nt-1)],[],1);
y = reshape(y+[s:s:s*nt],[],1);

% Detection of the unneeded non zero elements:
ind = x<=s*nt & y<=s*nt;

% Sparse matrix creation:
S = sparse(x(ind),y(ind),1,s*nt,s*nt);

% Plot the results:
spy(S)

这里我们只计算非零值的位置。我们避免预分配会减慢计算速度的大矩阵。

基准：

我已经使用 matlab 在线运行基准测试，可用内存有限。如果有人运行在他的本地计算机上进行具有更大价值的基准测试，请随意这样做。

通过这些配置，似乎使用隐式扩展确实更快。

基准代码：

for ii = 1:100
    s   = ii; k = 4; nt = ii;
    Auw = rand(s,s*k)>0.75;

    f_expa = @() func_expansion(s,nt,Auw);
    f_loop = @() func_loop(s,k,nt,Auw);

    t_expa(ii) = timeit(f_expa);
    t_loop(ii) = timeit(f_loop);
end

plot(1:100,t_expa,1:100,t_loop)
legend('Implicit expansion','For loop')
ylabel('Runtime (s)')
xlabel('x and nt value')

% obchardon suggestion
function S = func_expansion(s,nt,Auw)
    [x,y] = find(Auw);
    x = reshape(x+[0:s:s*(nt-1)],[],1);
    y = reshape(y+[s:s:s*nt],[],1);
    ind = x<=s*nt & y<=s*nt;
    S = sparse(x(ind),y(ind),1,s*nt,s*nt);
end

% Dev-il suggestion
function B = func_loop(s,k,nt,Auw)
    B = false(nt*s, (nt+k)*s);
    for it = 1:nt
        B( (it-1)*s+1:it*s, it*s+1:(it+k)*s ) = Auw;
    end
    B = B(1:nt*s, 1:nt*s);
end

通过 3D 数组页面的交错复制构建邻接矩阵

Construction of an adjacency matrix through staggered replication of pages of a 3D array

performance

matlab

matrix

graph-theory

adjacency-matrix

背景

设置

备注

演示 + 简单的实现