C 中克罗内克积的高效计算

Question

我对 C 语言相当陌生，在我的大部分研究中都不需要比 python 更快的东西。然而，事实证明，我最近所做的工作需要计算相当大的向量/矩阵，因此可能需要 C+MPI 解决方案。

从数学上讲，这个任务很简单。我有很多维数约为 40k 的向量，希望计算这些向量的选定对的 Kronecker Product，然后对这些 kronecker 乘积求和。

问题是，如何有效地做到这一点？下面的代码结构，使用for循环，或者获取效果有什么问题吗？

下面描述的函数kron 传递长度为vector_size 的向量A 和B，并计算它们的kronecker 乘积，并将其存储在C 中，即vector_size*vector_size 矩阵。

void kron(int *A, int *B, int *C, int vector_size) {

    int i,j;

    for(i = 0; i < vector_size; i++) {
        for (j = 0; j < vector_size; j++) {
            C[i*vector_size+j] = A[i] * B[j];
        }
    }
    return;
}

这对我来说似乎很好，当然（如果我没有犯一些愚蠢的语法错误）会产生正确的结果，但我偷偷怀疑嵌入的 for 循环不是最佳的。如果还有其他方法我应该解决这个问题，请告诉我。欢迎提出建议。

感谢您的耐心和任何建议。再说一次，我对 C 语言非常缺乏经验，但谷歌搜索并没有给我带来这个查询的乐趣。

Answer 1

为了解决你的问题，我认为你应该尝试使用 Eigen 3，它是一个使用所有矩阵函数的 C++ 库！

如果你有时间，去看看它的文档吧！ =)

祝你好运！

Answer 2

另一个易于实现的优化是，如果您知道数组的内部维度可以被 n 整除，则将 n 个赋值语句添加到循环体中，从而减少必要的迭代次数，并进行相应的更改循环计数。

可以通过在外部循环周围使用 switch 语句来概括这种策略，其中数组大小可以被 2、3、4 和 5 整除，或者最常见的任何情况。这可以带来相当大的性能提升，并且与建议 1 和 3 兼容，以进一步优化/并行化。一个好的编译器甚至可以为你做这样的事情（也就是循环展开）。

另一个优化是使用指针算法来避免数组索引。像这样的东西应该可以解决问题：

int i, j;

for(i = 0; i < vector_size; i++) {
    int d = *A++;
    int *e = B;

    for (j = 0; j < vector_size; j++) {
        *C++ = *e++ * d;
    }
}

这也避免了多次访问 A[i] 的值，方法是将其缓存在一个局部变量中，这可能会给您带来轻微的速度提升。 （请注意，此版本不可并行化，因为它会更改指针的值，但仍可用于循环展开。）

Answer 3

由于您的循环体都是完全独立的，因此肯定有一种方法可以加速这一过程。在考虑 MPI 之前，最简单的方法是利用多个内核。 OpenMP 在这方面应该做得很好。

#pragma omp parallel for
for(int i = 0; i < vector_size; i++) {
    for (int j = 0; j < vector_size; j++) {
        C[i][j] = A[i] * B[j];
    }
}

现在很多编译器都支持这一点。

您也可以尝试将一些常用表达式拖出内部循环，但体面的编译器，例如 gcc、icc 或 clang 应该可以自己完成此操作：

#pragma omp parallel for
for(int i = 0; i < vector_size; ++i) {
    int const x = A[i];
    int * vec = &C[i][0];
    for (int j = 0; j < vector_size; ++j) {
        vec[j] = x * B[j];
    }
}

顺便说一句，使用int 进行索引通常不是正确的做法。 size_t 是正确的 typedef 用于与对象的索引和大小有关的所有内容。

Answer 4

找到解决方案（感谢@Jeremiah Willcock）：GSL's BLAS bindings 似乎做得很好。如果我们逐步选择向量对 A 和 B 并将它们添加到一些“运行总”向量/矩阵 C，则上述 kron 函数的以下修改版本

void kronadd(int *A, int *B, int *C, int vector_size, int alpha) {

    int i,j;

    for(i = 0; i < vector_size; i++) {
        for (j = 0; j < vector_size; j++) {
            C[i*vector_size+j] = alpha * A[i] * B[j];
        }
    }
    return;
}

从功能上讲，精确对应于 BLAS DGER 功能（可通过gsl_blas_dger 访问）。初始的 kron 函数是 DGER，其中 alpha = 0 和 C 是正确维度的未初始化（归零）矩阵/向量。

事实证明，最终，对这些库简单地使用 python 绑定可能会更容易。但是，我想我在尝试解决这些问题时学到了很多东西。其他回复中有一些更有用的建议，如果您有同样的问题需要处理，请查看它们。谢谢大家！

Answer 5

如果您的编译器支持 C99（并且您从不传递与 A 和 B 相同的向量），请考虑在支持 C99 的模式下编译并将您的函数签名更改为：

void kron(int * restrict A, int * restrict B, int * restrict C, int vector_size);

restrict 关键字向编译器保证A、B 和C 指向的数组不会别名（重叠）。使用您编写的代码，编译器必须在每次执行内部循环时重新加载A[i]，因为它必须是保守的并假设您对C[] 的存储可以修改A[] 中的值。在restrict 下，编译器可以假设这不会发生。

Answer 6

这是数值计算界的一个常见问题，最好的办法是使用经过良好调试的包，如Matlab（或其Free Software clones之一）。

你甚至可以找到一个python binding，这样你就可以摆脱 C。

以上所有内容（可能）都会比严格用 python 编写的代码更快。如果您需要比这更快的速度，我建议您做几件事：

1. 考虑使用 Fortran 而不是 C。Fortran 编译器往往更擅长优化数值计算（如果您使用 gcc，则有一个例外，因为它的 C 和 Fortran 编译器使用相同的后端）。
2. 考虑并行化您的算法。我知道有一些 Fortran 变体具有并行循环语句。我认为周围有一些 C 插件可以做同样的事情。如果您使用的是 PC（和单精度），您还可以考虑使用视频卡的 GPU，它本质上是一个非常便宜的阵列处理器。

Answer 7

对于双精度向量（单精度和复数相似），您可以使用 BLAS 例程 DGER（排名一更新）或类似的方法一次做一个产品，因为它们是全部在向量上。你要乘多少个向量？请记住，添加一堆向量外积（您可以将 Kronecker 积视为）最终会作为矩阵-矩阵乘法，BLAS 的DGEMM 可以有效地处理。不过，如果您确实需要整数运算，您可能需要编写自己的例程。