Eigen C++中的按列点积

Question

有没有一种简单的方法来计算两个矩阵（称为A和B，类型为Eigen::MatrixXd）的列向量点积，这些矩阵的维度为mxn，而不必评估A*B或不必使用for循环？结果向量需要具有1xn或nx1的维度。此外，我正在尝试使用C++中的Eigen进行此操作。

- Zedd

---

逐元素相乘，然后求和？在MATLAB中，它将是sum(A .* B)。Eigen提供了这些操作，但我不知道调用的确切名称。 - Ben Voigt

不错！那应该可以了。谢谢！ - Zedd

可以使用“Eigen :: Map”将矩阵重塑为向量，然后取它们的内积。 - vsoftco

@Zedd：如果我的评论为您提供了制作工作代码所需的提示，请回答您自己的问题，展示使用Eigen函数完成它的方法，以便未来的访问者可以参考。 - Ben Voigt

Answer 1

有很多方法可以实现这个，它们都执行惰性求值：

res = (A.array() * B.array()).colwise().sum();
res = (A.cwiseProduct(B)).colwise().sum();

还有我的最爱：

res = (A.transpose() * B).diagonal();

- ggael

---

一个好奇的问题：通过将对角线附加到点积操作中，Eigen是否只会在内部计算（A的第i行）与（B的第i列）的点积，以避免不必要的计算？对于这三种实现，哪一种会更快？（因为第一种和第二种都会避免不必要的计算）。 - JackieLam

正如答案中所写，它们都执行惰性求值，因此最后一个仅计算对角线条目。它们应该生成大致相同的代码。您可以在编译器资源管理器上进行检查。 - ggael

(A.transpose() * B).diagonal(); 这个表达式的渐近性能是否正确？同时，为什么 A.colwise().dot(B.colwise()) 不起作用？ - Alec Jacobson

Answer 2

我基于@ggael的答案进行了实验。

MatrixXd A = MatrixXd::Random(600000,30);
MatrixXd B = MatrixXd::Random(600000,30);

MatrixXd res;
clock_t start, end;
start = clock();
res.noalias() = (A * B.transpose()).diagonal();
end = clock();
cout << "Dur 1 : " << (end - start) / (double)CLOCKS_PER_SEC << endl;

MatrixXd res2;
start = clock();
res2 = (A.array() * B.array()).rowwise().sum();
end = clock();
cout << "Dur 2 : " << (end - start) / (double)CLOCKS_PER_SEC << endl;

MatrixXd res3;
start = clock();
res3 = (A.cwiseProduct(B)).rowwise().sum();
end = clock();
cout << "Dur 3 : " << (end - start) / (double)CLOCKS_PER_SEC << endl;

输出结果为：

Dur 1 : 10.442
Dur 2 : 8.415
Dur 3 : 7.576

看起来，diagonal()解决方案是最慢的。cwiseProduct是最快的。 而且内存使用相同。

- JackieLam

---

@ggael，请随意发表评论。 - JackieLam

我进行了一个非常类似的测试，但使用了16x3个单精度矩阵，并且每个操作都执行了一百万次。使用 clang++ -O3，我得到了175、129、131毫秒的时间，这与您的结果相似。 - wcochran

Answer 3

这是我使用 Eigen::Map 的方法（假设矩阵为实数，可通过取伴随转换为复数），其中 rows 和 cols 表示行数和列数：

#include <Eigen/Dense>
#include <iostream>

int main()
{
    Eigen::MatrixXd A(2, 2);
    Eigen::MatrixXd B(2, 2);
    A << 1, 2, 3, 4;
    B << 5, 6, 7, 8;

    int rows = 2, cols = 2;

    Eigen::VectorXd vA = Eigen::Map<Eigen::VectorXd>(
                             const_cast<double *>(A.data()), rows * cols, 1);
    Eigen::VectorXd vB = Eigen::Map<Eigen::VectorXd>(
                             const_cast<double *>(B.data()), rows * cols, 1);

    double inner_prod = (vA.transpose() * vB).sum();

    std::cout << inner_prod << std::endl;
}

- vsoftco

---

我现在理解您想要一个向量作为结果？也就是说，您需要一个\vect{result}，其中每个分量都是A(i,)*B(:,i)'吗？ - vsoftco

没错。我不确定这如何转换成你上面的代码。 - Zedd