在C语言中高效地计算Kronecker积

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语言非常不熟悉，但在谷歌上搜索并没有为我的问题带来太多的帮助。

- Edward Grefenstette

---

两个向量的 Kronecker 乘积不是一个向量吗？ - BlueRaja - Danny Pflughoeft

一个好问题加1。欢迎来到SO。 - Jens Gustedt

Jens: 谢谢 :-) Dan: 一些库（如 scipy，我认为还有 matlab）会将两个 N 维向量的外积得到的 N x N 矩阵展开成一个 N*N 维向量。我不太在意我得到的是矩阵还是向量形式，只要我能对其进行求和即可... - Edward Grefenstette

可能是一个小错误：按照现在的写法，C[i][j] 无法编译通过。你需要像这样 C[i*vector_size+j] 或者声明函数参数为 int C[vector_size][]（使用gcc重新排列）。 - Joseph Quinsey

@Joseph：不是打错字，而是更多证明我通常是一个Python的“码农” ：P 感谢你的纠正！ - Edward Grefenstette

我也犯了一个错别字：第二个建议应该是 C[][vector_size]。 - Joseph Quinsey

Answer 1

由于您的循环体完全独立，因此肯定有一种加速方法。在考虑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。

- Jens Gustedt

---

这听起来很有趣。我会稍微研究一下openMP。（顺便说一句：我不是在谈论针对这个特定任务使用MPI，而是更多地针对并行计算不同向量的问题——这是一个单独的问题）。冒昧问一句，您能否详细解释一下“您还可以尝试从内部循环中拖出一些常见表达式”的意思？ - Edward Grefenstette

谢谢！我不知道那些技巧。非常有帮助 :-) - Edward Grefenstette

Answer 2

对于双精度向量（单精度和复数类似），您可以使用BLAS例程DGER（秩一更新）或类似方法逐个执行乘积，因为它们都是向量。 您要乘多少个向量？请记住，添加一堆向量外积（您可以将Kronecker积视为此）最终变成矩阵-矩阵乘法，BLAS的DGEMM可以高效地处理。 如果您确实需要整数操作，则可能需要编写自己的例程。

- Jeremiah Willcock

---

BLAS是我在探索中遇到的东西。然而，我在我的实验室机器上使用它时遇到了很多麻烦（似乎无法将cblas.h放入Fedora Core），甚至连基本教程都无法完成。我很难找到可以理解的文档。我认为使用它和记录它的人比我操作的水平稍高 :-P - Edward Grefenstette

@egrefen：GSL（http://www.gnu.org/software/gsl/）可能是一个易于安装的软件包；还有 Goto BLAS （http://www.tacc.utexas.edu/tacc-projects/gotoblas2/）和 ATLAS（http://math-atlas.sourceforge.net/）。如果需要供应商定制版本，则AMD有他们的ACML，英特尔则有MKL。 - Jeremiah Willcock

谢谢，这很有帮助。那里还有更多的文档。我会看一下它是否符合我的需求... - Edward Grefenstette

Answer 3

如果您的编译器支持C99（且您从未将相同的向量作为A和B传递），请考虑以支持C99的模式进行编译，并将函数签名更改为：
A和B不再是相同的向量。

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

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

- caf

Answer 4

解决方案找到了（感谢@Jeremiah Willcock）：GSL的BLAS绑定似乎非常适合这个问题。如果我们正在逐步选择向量对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函数是带有alpha = 0和C为正确维度的未初始化（清零）矩阵/向量的DGER。

事实证明，最终使用这些库的Python绑定可能会更容易。不过，在尝试弄清楚这些东西的过程中，我觉得我学到了很多。如果您遇到类似的问题，请查看其他回复中的一些有用建议。感谢大家！

- Edward Grefenstette

Answer 5

另一个容易实现的优化是，如果您知道数组的内部维度可以被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]的值多次访问它，这可能会给您带来轻微的速度提升。（请注意，此版本不可并行化，因为它改变了指针的值，但仍可以使用循环展开。）

- Keith

Answer 6

在数值计算方面，这是一个常见的问题，最好的方法是使用像 Matlab（或其 自由软件克隆版 之一）这样的经过良好调试的包。

你甚至可以找到一个 Python绑定 版本，这样就可以摆脱C语言了。

以上所有方法（可能）都比纯Python编写的代码更快。如果你需要更高的速度，我建议采取以下几种方法：

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

- T.E.D.

---

是的，Python有很多好用的库（如numpy、scipy），可以轻松处理这种情况。然而，它们并不像C语言那样高效，也不太便携（比如无法在我们大学的超级计算机设施上使用）。另外，我想亲自动手用C语言进行学习。与已经拥有高效实现的Matlab相比，C语言更具吸引力的一点是，它可以成为一个轻量级的Python扩展，并且我们可以在GPL或类似的许可下发布最终的框架。感谢你的建议。 - Edward Grefenstette

@Edward Grefenstette - 好的...我的回复有点长，所以我把它加到了我的答案中。 - T.E.D.

Answer 7

uint32_t rA  = 3;
uint32_t cA  = 5;
uint32_t lda = cA;
uint32_t rB  = 5;
uint32_t cB  = 3;
uint32_t ldb = cB;
uint32_t rC  = rA*rB;
uint32_t cC  = cA*cB;
uint32_t ldc = cC;
double *A = (double *)malloc(rA*cA*sizeof(double));
double *B = (double *)malloc(rB*cB*sizeof(double));
double *C = (double *)malloc(rC*cC*sizeof(double));
for (uint32_t i=0, allA=rA*cA; i<allA; i++)
    A[i]=i;
for (uint32_t i=0, allB=rB*cB; i<allB; i++)
    B[i]=i;
for (uint32_t i=0, allC=rC*cC; i<allC; i++)
    C[i]=0;
for (uint32_t i=0, allA=rA*cA; i<allA; i++)
{
    for (uint32_t j=0, allB=rB*cB; j<allB; j++)
      C[((i/lda)*rB+j/ldb)*ldc
       + (i%lda)*cB+j%ldb     ]=A[i]*B[j];
}

- OChicken

Answer 8

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

如果有时间，去看看它的文档！=）

祝你好运！

- Sabrina Tesla