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

Question

背景

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

环境

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

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

第三维的含义是转换需要多少时间步，例如：如果A(1,3,2)非零，则表示状态#1可以转换为状态#3，这将需要2时间步。

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

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

其中主块对角线由nt0 个块组成，并且 的切片A逐渐“向右推”直到“时间用完”， 的切片A最终“超出”了B? 表明不可能有更多的转换。由于B由nt*nt s*s块组成，因此其大小为(nt*s)×(nt*s)。

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

笔记

• 由于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在线运行基准测试，可用内存有限。如果有人想在本地计算机上运行具有更大价值的基准测试，请随意这样做。

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

基准代码：

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