时间戳

首页 > 解决方案 > 通过 3D 数组页面的交错复制构建邻接矩阵

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

问题描述

背景

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

环境

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

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]"

其中主块对角线由nt0 块组成，并且 切片A逐渐“推”到右侧直到“时间用完”，最终切片A“超出” B⇒ 表示不再可能进行转换。由于B由nt*nt s*s块组成，其大小为(nt*s)×(nt*s).

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

笔记

演示 + 朴素的实施

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);

导致：

在此处输入图像描述

标签： performancematlabmatrixgraph-theoryadjacency-matrix

解决方案

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

% 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

推荐阅读