计算矩阵每个元素的指数的最有效方法答案

【问题标题】：Most efficient way to calculate the exponential of each element of a matrix计算矩阵每个元素的指数的最有效方法
【发布时间】：2010-07-23 21:31:42
【问题描述】：

我正在从 Matlab 迁移到 C + GSL，我想知道计算矩阵 B 的最有效方法是什么：

B[i][j] = exp(A[i][j])

其中 i 在 [0, Ny] 中，j 在 [0, Nx] 中。

请注意，这与矩阵指数不同：

B = exp(A)

这可以通过 GSL (linalg.h) 中的一些不稳定/不受支持的代码来实现。

我刚刚找到了蛮力解决方案（几个“for”循环），但是有什么更聪明的方法吗？

编辑

Drew Hall 的解答帖结果

所有结果都来自一个 1024x1024 for(for) 循环，其中在每次迭代中分配了两个 double 值（一个复数）。 时间是 100 次执行的平均时间。

考虑到存储矩阵的 {Row,Column}-Major 模式时的结果：
- 在 Row-Major 模式下循环内循环中的行时为 226.56 毫秒（案例 1）。
- 在 Row-Major 模式下在内循环中循环列时为 223.22 毫秒（案例 2）。
- 使用 GSL 提供的gsl_matrix_complex_set 函数时为 224.60 毫秒（案例 3）。

案例1的源代码：

for(i=0; i<Nx; i++)
{
    for(j=0; j<Ny; j++)
    {
        /* Operations to obtain c_value (including exponentiation) */
        matrix[2*(i*s_tda + j)] = GSL_REAL(c_value);
        matrix[2*(i*s_tda + j)+1] = GSL_IMAG(c_value);
    }
}

案例2的源代码：

for(i=0; i<Nx; i++)
{
    for(j=0; j<Ny; j++)
    {
        /* Operations to obtain c_value (including exponentiation) */
        matrix->data[2*(j*s_tda + i)] = GSL_REAL(c_value);
        matrix->data[2*(j*s_tda + i)+1] = GSL_IMAG(c_value);
    }
}

案例3的源代码：

for(i=0; i<Nx; i++)
{
    for(j=0; j<Ny; j++)
    {
        /* Operations to obtain c_value (including exponentiation) */
        gsl_matrix_complex_set(matrix, i, j, c_value);
    }
}

【问题讨论】：

我认为您只是认为它效率低下/愚蠢，因为您已经习惯了 MATLAB，其中使用嵌套 fors 是灾难的根源。但是 MATLAB 会在幕后做到这一点，只是如果你把它写出来，那就糟透了。但是 C 不是那样工作的，你总是把它写出来，它会和 MATLAB 一样快/快...
嗯...这些性能数据无法区分。我很惊讶——我们能看到你在每种情况下运行的代码吗？我想知道您的编译器是否注意到这些值是独立的并自动重新排序和展开循环以最大化性能。在每种情况下查看生成的程序集都会告诉我们这一点。
在每种情况下，您都在赋值的 RHS 上使用了一个常量值（exp() 发生了什么？）——我怀疑编译器已将其提升到循环之外并替换了循环使用某种类似 memset() 的操作。尝试根据数据进行分配——即使只是像 matrix[2*(is_tda + j)] = sin(matrix[2*(is_tda + j)]) 之类的东西也可能击败优化并显示缓存效果。
@Drew：c_value 不是一个常数。我只是缩短了代码以避免它太长。我使用的实际代码托管在signal.c第45行附近的项目gico-lib（有点混乱）的google代码中。我将尝试使用该技巧伪造编译器（达尔文的 gcc 4.2.1 (MacOsX)）并稍后发布结果。
涉及到矩阵元素的赋值（A[i] = exp(A[i]) 等价）情况1和2的差别差不多，大约5毫秒。

标签： c performance matrix exponential

【解决方案1】：

没有办法避免遍历所有元素并在每个元素上调用exp() 或等效项。但是有更快和更慢的迭代方式。

特别是，您的目标应该是尽量减少缓存未命中。找出您的数据是以行优先还是列优先的顺序存储的，并确保安排循环，以便内部循环遍历内存中连续存储的元素，而外部循环则大步前进到下一行（如果行主要）或列（如果列主要）。虽然这看起来微不足道，但它可以在性能上产生巨大的差异（取决于矩阵的大小）。

处理完缓存后，您的下一个目标是消除循环开销。第一步（如果您的矩阵 API 支持它）是从嵌套循环（M 和 N 边界）转到遍历基础数据的单个循环（MN 边界）。为此，您需要获取指向底层内存块的原始指针（即双精度而不是双精度**）。

最后，加入一些循环展开（即，为循环的每次迭代执行 8 或 16 个元素）以进一步减少循环开销，这可能会尽可能快。您可能需要一个带有贯穿的最终 switch 语句来清理剩余元素（当您的数组大小 % 块大小！= 0 时）。

【讨论】：

哇，这是我之前没有想到的策略变化。谢谢，我会把它写下来并尝试将它应用到我所做的许多嵌套 fors 循环中。
当然。完成后，请在此处添加另一条评论（或编辑您的问题），让我们知道您是否确实获得了任何性能改进！
他将如何确定矩阵是以行优先还是列顺序存储？我希望自己能够使用这种优化。另外，Drew Hall = for loop superhero
@Rafe：通常具有 C 遗产的代码是行主要的，而具有 Fortran 遗产的代码（或与 LAPACK 等 Fortran 库的接口）是列主要的。查看您正在使用的矩阵库的文档，或者如果所有其他方法都失败了，请使用非常大的矩阵和时间尝试两种方式 - 它应该很明显。
我编辑了问题以包含我获得的一些结果。稍后我会尝试其他优化。

【解决方案2】：

不，除非有一些我没有听说过的奇怪的数学怪癖，否则您几乎只需要使用两个 for 循环遍历元素即可。

【讨论】：

身边有魔法数学家吗？做双 for 循环确实看起来效率低下。
Alejandro：如果你想要N_x * N_y 独立计算指数，那么A 不再是矩阵，而只是一堆数字。 :-)
如果你想要数学家，我会去数学论坛或类似性质的东西（如 www.askdrmath.com）并尝试他们的知识。试一试也无妨！
kaize.se：你说得对，闻起来不像改进双 for 循环 :-( 顺便说一句，我按照 Rafe 的建议将问题发送给了 Math 博士，以防万一某个神奇的数学家在线。

【解决方案3】：

如果您只想将exp 应用于数字数组，那么真的没有捷径可走。你必须调用它 (Nx * Ny) 次。如果一些矩阵元素很简单，比如 0，或者有重复的元素，一些记忆可能会有所帮助。

但是，如果您真正想要的是矩阵指数（这非常有用），我们依赖的算法是DGPADM。它在 Fortran 中，但您可以使用 f2c 将其转换为 C。Here's the paper on it.

【讨论】：

【解决方案4】：

由于没有显示循环的内容，计算c_value的位我们不知道代码的性能是受内存带宽限制还是受CPU限制。唯一确定的方法是使用分析器，并且是一个复杂的分析器。它需要能够测量内存延迟，即 CPU 空闲等待数据从 RAM 到达的时间量。

如果您受到内存带宽的限制，那么一旦您按顺序访问内存，您就无能为力了。当按顺序获取数据时，CPU 和内存工作得最好。随机访问会影响吞吐量，因为数据更有可能必须从 RAM 中提取到缓存中。您总是可以尝试获得更快的 RAM。

如果您受到 CPU 的限制，那么您可以使用更多选项。使用 SIMD 是一种选择，手动编码浮点代码也是如此（由于许多原因，C/C++ 编译器在 FPU 代码方面表现不佳）。如果这是我，并且内部循环中的代码允许这样做，我将有两个指向数组的指针，一个在开始处，另一个在 4/5 处通过它。每次迭代，将使用第一个指针执行 SIMD 操作，使用第二个指针执行标量 FPU 操作，以便循环的每次迭代执行五个值。然后，我将 SIMD 指令与 FPU 指令交错以降低延迟成本。这应该不会影响您的缓存，因为（至少在 Pentium 上）MMU 可以同时传输多达四个数据流（即在没有任何提示或特殊说明的情况下为您预取数据）。

【讨论】：