eigen 是否有像 H.transpose()*H 这样的自转置乘法优化答案

【问题标题】：does eigen have self transpose multiply optimization like H.transpose()*Heigen 是否有像 H.transpose()*H 这样的自转置乘法优化
【发布时间】：2023-03-31 05:56:01
【问题描述】：

我浏览了 eigen 的教程 https://eigen.tuxfamily.org/dox-devel/group__TutorialMatrixArithmetic.html

它说 “注意：对于担心性能的 BLAS 用户，诸如 c.noalias() -= 2 * a.adjoint() * b; 之类的表达式已完全优化并触发单个类似 gemm 的函数调用。”

但是像 H.transpose() * H 这样的计算怎么样，因为它的结果是一个对称矩阵，所以它应该只需要正常 A*B 的一半时间，但在我的测试中，H.transpose() * H 花费相同time as H.transpose() * B. eigen 是否对这种情况有特殊的优化，和opencv一样，有类似的功能。

我知道对称优化会破坏向量化，我只是想知道本征是否有解决方案可以同时提供对称优化和向量化

【问题讨论】：

标签： optimization eigen matrix-multiplication neon

【解决方案1】：

你是对的，你需要告诉 Eigen 结果是这样对称的：

Eigen::MatrixXd H = Eigen::MatrixXd::Random(m,n);
Eigen::MatrixXd Z = Eigen::MatrixXd::Zero(n,n);
Z.template selfadjointView<Eigen::Lower>().rankUpdate(H.transpose());

最后一行在下三角部分计算Z += H * H^T。上半部分保持不变。你想要一个完整的矩阵，然后将下部复制到上部：

Z.template triangularView<Eigen::Upper>() = Z.transpose();

这个rankUpdate 例程是完全矢量化的，可与等效的 BLAS 相媲美。对于小矩阵，更好地执行完整的产品。

另请参阅各自的doc。

【讨论】：

这绝对是我想要的，我已经测试它正确，但在我的测试中我认为 Z.sefladjointView().rankUpdate(H);意思是 Z+= H*H' 我是对的吗？
对，如果你想要H'*H，那就打电话给.rankUpdate(H.adjoint());。
很抱歉再次打扰您:) 我在计算 HAH' 时遇到了问题，其中 A = A' ，那么有什么方法可以加速它吗？跨度>
你可以像往常一样计算B=H*A，然后Z.triangularView<Lower>() = B*H'。这将只计算第二个产品的一半。
有没有考虑使用.rankUpdate(H.adjoint()); vs .rankUpdate(H.transpose());？从经验试验来看，它们都产生了相同的结果，并且使用 H 的随机数据似乎具有非常相似的速度性能。