performance - 通过交错复制 3D 阵列的页面构建邻接矩阵

Question

背景
我正在尝试建模一个可以在每个时间步更改其配置的系统。各种配置是预先知道的，不依赖于时间步长。在某些配置之间允许转换，而在其他配置之间禁止转换。目标是构建一个跨越多个时间步长的允许转换的邻接矩阵。
环境
让 A成为 s*s*k表示允许转换的逻辑矩阵，以及 A1...Ak代表 A 的页面/切片:

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

第三维的含义是一个过渡需要多少个时间步，例如:if 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 的切片逐渐向右“推”直到“时间用完”， A的切片最终在 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 .

欢迎使用外部脚本/类/工具。

我的一个想法不是构建最初交错的矩阵，而是具有 [A1] 的主对角线块和 circshift -ing 和掩蔽，当其他一切都完成时。

演示 + Naïve 实现

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运行基准测试，可用内存有限。如果有人会在他的本地计算机上以更大的值(value)运行基准测试，请随意这样做。

使用这些配置，使用隐式扩展似乎确实更快。
基准代码:

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