> 4D阵列MatLab的巨大性能损失原因？

Question

简介

我有一个算法可以循环数十亿（万亿）次并操作存储在 7 维 [10x10x10x10x10x10x10] 中的矩阵，我发现访问 7 维矩阵中的元素非常慢，我很好奇我运行了一些测试确定访问多维矩阵元素的性能。

假设

有人提醒我，MatLab 在后台使用线性索引，我的一个朋友说性能损失可能是由于将“正常”索引转换为后台Source 的线性索引。

测试方法

为了验证这个假设，我测试了使用线性索引和常规索引访问 2D 到 7D 矩阵的元素。我改变了我正在访问的元素以及我正在访问的矩阵大小，即每个维度的长度，但这并没有显着改变结果。我用于测试的文件在下面找到。使用的硬件是 Intel(R) Xeon(R) CPU E5-1620 v2 @ 3.70 GHz，16GB RAM。 MatLab 版本为 R2013B。

结果

Normal indexing 2D: 0.073659s
Linear indexing 2D: 0.064026s
Normal indexing 3D: 0.050719s
Linear indexing 3D: 0.064096s
Normal indexing 4D: 0.055674s
Linear indexing 4D: 0.062112s
Normal indexing 5D: 15.689907s
Linear indexing 5D: 5.265076s
Normal indexing 6D: 16.660496s
Linear indexing 6D: 5.295958s
Normal indexing 7D: 17.029072s
Linear indexing 7D: 5.291139s

在处理 4D 矩阵的性能方面，与线性索引相比，普通索引似乎非常相似。对于 3D 和 4D 矩阵，使用普通索引似乎稍微好一些。在 4D 矩阵之上，线性索引比正常索引更可取，但正常索引和常规索引都会受到巨大的性能损失（约两个数量级）。

结论

这个故事的寓意是在以高性能为目标时，在运行 MatLab 时仔细考虑矩阵中需要超过四个维度（除了显而易见的事实，例如 C++ 对于大多数应用程序来说要快得多等等，但这是另一个讨论可能）。

问题

如标题所示：当矩阵中超过四个维度时，性能大幅下降的（潜在）原因是什么？其他语言（例如 C++）是否也表现出这种行为？

用于测试的代码

clear all
clc

% Number if iterations 
n=10000000;

A = rand(10,10);
k1=sub2ind(size(A),3,2);
fprintf('\n')
tic
for ii=1:n
    A1=A(3,2);
end
a=toc;
fprintf('Normal indexing 2D: %fs\n',a)


tic
for ii=1:n
    A2=A(k1);
end
a=toc;
fprintf('Linear indexing 2D: %fs\n',a)

B = rand(10,10,10);
k2=sub2ind(size(B),3,2,1);

tic
for ii=1:n
    B1=B(3,2,1);
end
a=toc;
fprintf('Normal indexing 3D: %fs\n',a)
tic
for ii=1:n
    B2=B(k2);
end
a=toc;
fprintf('Linear indexing 3D: %fs\n',a)

C = rand(10,10,10,10);
k3=sub2ind(size(C),3,2,1,10);

tic
for ii=1:n
    C1=C(3,2,1,10);
end
a=toc;
fprintf('Normal indexing 4D: %fs\n',a)
tic
for ii=1:n
    C2=C(k3);
end
a=toc;
fprintf('Linear indexing 4D: %fs\n',a)

D = rand(10,10,10,10,10);
k4=sub2ind(size(D),3,2,1,10,1);

tic
for ii=1:n
    D1=D(3,2,1,10,1);
end
a=toc;
fprintf('Normal indexing 5D: %fs\n',a)
tic
for ii=1:n
    D2=D(k4);
end
a=toc;
fprintf('Linear indexing 5D: %fs\n',a)

E = rand(10,10,10,10,10,10);
k5=sub2ind(size(E),3,2,1,10,1,2);

tic
for ii=1:n
    E1=E(3,2,1,10,1,2);
end
a=toc;
fprintf('Normal indexing 6D: %fs\n',a)
tic
for ii=1:n
    E2=E(k5);
end
a=toc;
fprintf('Linear indexing 6D: %fs\n',a)

F = rand(10,10,10,10,10,10,10);
k6=sub2ind(size(F),3,2,1,10,1,2,3);

tic
for ii=1:n
    F1=F(3,2,1,10,1,2);
end
a=toc;
fprintf('Normal indexing 7D: %fs\n',a)
tic
for ii=1:n
    F2=F(k6);
end
a=toc;
fprintf('Linear indexing 7D: %fs\n',a)

Answer 1

在某些时候，一个值必须有一个线性地址，因为这是大多数处理器的工作方式。处理器非常擅长从单个位置获取数据。

在整个答案中，术语 matrix 和 array 将互换使用。

单维矩阵
一维数组的地址计算：

  address = starting_array_address + (index * sizeof(object_type));

object_type是对象的类型，如unsigned char或int。

二维矩阵
以一维数组表示的二维数组中某个值的索引位置的计算为：

  index = Dimension2_value * Dimension1_capacity + Dimension1_value;

值的地址可以使用一维对象的索引和方程来计算。

3 维矩阵
Dimension3_value 可以表示为立方体内的平面数。
因此，这将等式扩展了另一个术语：

index = Dimension_3_value * (Dimension2_capacity)
      + Dimension2_value * Dimension1_capacity
      + Dimension1_value;

4D 及以上矩阵 4D 及以上矩阵的计算可以以类似的方式进行扩展。

索引计算的瓶颈是乘法。

某些平台可以并行计算项，因此具有一定的速度优势。

速度推测
我猜 Matlab 已经优化了 2D 和 3D 矩阵的索引计算，但没有优化更大的矩阵。 2D 和 3D 矩阵比 4D 和更大的矩阵使用更广泛。

典型的优化是将每个项的值分配给不同的寄存器，然后对寄存器求和。这是为了使值尽可能接近处理器。一些处理器可能没有足够的寄存器来进行 4D 计算，因此编译器可能不得不将值放在临时内存中（处理器外部），从而降低性能。

获取值
对于 1D、2D 和 3D 矩阵，处理器可能能够将数据包含在其数据缓存中，从而提高性能。随着维度的增长，矩阵的大小必须缩小以适应处理器的数据缓存。任何时候一个值不在缓存中都称为缓存未命中，处理器必须从处理器外部的内存中加载该值。根据重新加载数据的方式，它可能会重新加载数据“条”或更多。无论哪种情况，都会降低性能。

总结
对于矩阵的每个附加维度，计算索引需要更多的处理时间。乘法是计算的瓶颈（即使它们可能是一条指令，乘法比加法花费的时间更长）。虽然每一项指标计算都可以并行执行，但可能没有足够的线程或核心或处理单元来一次执行所有维度的指标计算。可能没有足够的寄存器或数据缓存来进行最佳计算，因此处理器将不得不从外部存储器中获取数据，从而降低性能。如果数据项不在处理器的数据高速缓存中，则处理器将浪费时间重新加载整个或部分数据高速缓存（因为处理器可能正在执行计算而不是重新加载高速缓存，所以浪费了时间）。

一般来说，对行列式（较小的子矩阵）进行操作通常会比在大型多维矩阵上进行操作获得更好的性能。