如果矩阵密集且不对称，特征库矩阵/向量操作是否比 .net 更快？答案

【问题标题】：Is Eigen library matrix/vector manipulation faster than .net ones if the matrix is dense and unsymmetrical?如果矩阵密集且不对称，特征库矩阵/向量操作是否比 .net 更快？
【发布时间】：2017-03-16 08:07:34
【问题描述】：

我有一些矩阵运算，主要处理诸如遍历矩阵的所有行和列并执行乘法运算a*mat[i,j]*mat[ii,j]：

public double[] MaxSumFunction()
{
   var maxSum= new double[vector.GetLength(1)];
   for (int j = 0; j < matrix.GetLength(1); j++)
   {
        for (int i = 0; i < matrix.GetLength(0); i++) 
        {
             for (int ii = 0; ii < matrix.GetLength(0); ii++)
              {
                   double wi= Math.Sqrt(vector[i]);
                   double wii= Math.Sqrt(vector[ii]);
                   maxSum[j] += SomePowerFunctions(wi, wii) * matrix[i, j]*matrix[ii, j];
              }
          }                      
     }
 }

    private double SomePowerFunctions(double wi, double wj)
    {

        var betaij = wi/ wj;
        var numerator = 8 * Math.Sqrt(wi* wj) * Math.Pow(betaij, 3.0 / 2)
            * (wi+ betaij * wj);
        var dominator = Math.Pow(1 - betaij * betaij, 2) +
            4 * wi* wj* betaij * (1 + Math.Pow(betaij, 2)) +
            4 * (wi* wi+ wj* wj) * Math.Pow(betaij, 2);


        if (wi== 0 && wj== 0)
        {
            if (Math.Abs(betaij - 1) < 1.0e-8)
                return 1;
            else
                return 0;
        }

        return numerator / dominator;
    }

如果矩阵很大，我发现这样的循环会特别慢。

我希望速度快。所以我正在考虑使用 Eigen 库重新实现这些算法。

我的矩阵不是对称的，不是稀疏的，并且不包含任何求解器都可以可靠利用的规律性。

我读到 Eigen solver can be fast 是因为：

编译器优化
矢量化
Multi-thread support

但我想知道考虑到我的矩阵特性，这些优势是否真的适用？

注意：我可以运行一两个样本来找出答案，但我相信在这里提出问题并将其记录在 Internet 上也会对其他人有所帮助。

【问题讨论】：

我不明白......如果你只需要性能，你为什么不使用一些低级别的东西，比如 OpenBLAS？为什么你需要重新发明人们花了几十年时间优化的矩阵乘法？您在那里标记的操作：mat[i,j]*mat[ii,j]，只是矩阵乘法与转置mat[j,ii]。任何 BLAS 接口库都可以为您做到这一点。顺便说一句，像这样循环是最慢的方法。如果您可以使用std::transform 执行此操作，您仍然可以获得更好的性能，这将为您启用矢量化。
@TheQuantumPhysicist, 1) mat[i,j]*mat[ii,j] 这不是转置操作 2)，我正在使用 C#，所以 std::transform 没有帮助 3) 我认为没有理由更喜欢 OpenBlas 而不是 Eigen，我很了解 Eigen，但 OpenBlas 不是这样。
我从没说过它是转置运算，我说它是与mat[ii,j]的转置相乘，所以它是mat[i,j]*Tr(mat[j,ii])，它意味着这个操作可以简化为矩阵乘法，并且由于矩阵乘法是本书中最古老的问题，你应该考虑使用 BLAS，因为它是线性代数的库。如果它是专门的 OpenBLAS（这是 BLAS 的一种实现）并不重要，但是，同样，你会在那里重新发明轮子。我遇到了大小超过 10000x10000 的矩阵乘法的量子力学问题，所以看看 BLAS。
@TheQuantumPhysicist，很好的解释。让我看看
@TheQuantumPhysicist，我不认为这可以直接转换为矩阵-矩阵乘积，因为求和不适用于公共索引j。另一方面，如果您将f(wi,wii) 视为一个矩阵，那么您可以做一些事情......此外，Eigen 提供了与其他优化的 BLAS 实现相似的性能量级，它还公开了一个 BLAS 兼容的 API，并且可以使用另一个 BLAS 实现....

标签： c++ .net matrix eigen

【解决方案1】：

在考虑低级优化之前，请查看您的代码并观察许多数量被多次重新计算。例如，f(wi,wii) 不依赖于j，因此它们可以预先计算一次（见下文），或者您可以重写循环以使 j 上的循环成为嵌套循环。然后嵌套循环将只是一个常数标量和矩阵的两列之间的系数乘积（我不是 .net 并假设 j 是索引列）。如果存储是列主要的，那么这个操作应该由你的编译器完全向量化（同样，我不知道.net，但任何 C++ 编译器都可以，如果你是 Eigen，它将被显式向量化）。这应该足以获得巨大的性能提升。

根据matrix 的大小，您还可以尝试通过将f(wi,wii) 预计算为MatrixXd F;（使用Eigen 的语言）来利用优化的矩阵-矩阵实现，然后观察整个计算量为：

VectorXd v = your_vector;
MatrixXd F = MatrixXd::nullaryExpr(n,n,[&](Index i,Index j) {
                 return SomePowerFunctions(sqrt(v(i)), sqrt(v(j)));
             });
MatrixXd M = your_matrix;
MatrixXd FM = F * M;
VectorXd maxSum = (M.array() * FM.array()).colwise().sum();

【讨论】：