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

Question

背景

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

设置

让A 是一个s*s*k 逻辑矩阵，表示允许的转换，A1...Ak 表示A 的页面/切片：

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

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

让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 的切片逐渐“推”到右侧直到“时间用完”，A 的切片最终“向外”的B ⇒ 表示不再可能进行转换。由于B 由nt*nt s*s 块组成，其大小为(nt*s)×(nt*s)。

问题：给定A 和nt，我们如何以最节省CPU 和内存的方式构造B？

注意事项

• 由于B 大部分都用零填充，所以sparse 可能是有意义的。
• 在我的应用程序中，CPU 效率（运行时）比内存效率更重要。
• 在真正的问题中，s=250 和 nt=6000。
• 欢迎使用外部脚本/类/工具。
• 我的一个想法不是最初构建交错矩阵，而是使用[A1] 块和circshift-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