如何在Eigen中优化稀疏与稠密自伴随矩阵的乘积运算

Eigen目前没有直接支持稀疏矩阵S的selfadjointView.rankUpdate(S,H)重载，但我们可以利用自伴随矩阵的特性，只计算结果的下三角部分并同步到J，从而节省约50%的计算量。

实现思路 #

1. 预计算稠密矩阵SH = S * H：稀疏矩阵乘稠密矩阵的结果为稠密矩阵，这一步是必要的，但后续计算可避免生成完整的SH*S^H矩阵。
2. 利用自伴随矩阵的对称性：(S*H*S^H)(i,j) = conj((S*H*S^H)(j,i))，因此只需计算i>=j的下三角元素，通过SelfAdjointView自动同步上三角，无需重复计算。
3. 通过向量点积高效计算下三角元素：利用Eigen的向量运算优化，避免手动遍历每个矩阵元素的低效操作。

具体代码 #

#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <complex>

using namespace Eigen;

Matrix<std::complex<double>, Dynamic, Dynamic> J, H;
SparseMatrix<std::complex<double>> S;

// ... 初始化H、J、S的代码（确保尺寸匹配，H、J为自伴随矩阵） ...

// 预计算S*H，得到稠密矩阵SH
MatrixXcd SH = S * H;
// 将S的伴随矩阵转为稠密形式，方便快速访问行数据
MatrixXcd S_adj_dense = S.adjoint();
// 获取J的下三角自伴随视图，修改下三角会自动同步上三角的共轭值
auto J_lower = J.selfadjointView<Lower>();

const int n = J.rows();
for (int i = 0; i < n; ++i) {
    for (int j = 0; j <= i; ++j) {
        // 计算(S*H*S^H)(i,j) = SH的第i行 与 S^H的第j行 的点积
        std::complex<double> elem = SH.row(i).dot(S_adj_dense.row(j));
        J_lower.coeffRef(i, j) += elem;
    }
}

性能优化与权衡 #

• 内存vs计算量：如果S的稀疏度极高，预存S_adj_dense会占用较多内存，此时可改为直接从稀疏矩阵提取列的共轭：VectorXcd s_col_conj = S.col(j).conjugate();，再与SH.row(i)做点积，以内存开销换取计算效率的平衡。
• 对比默认实现：直接使用J.noalias() += S*H*S.adjoint();会生成完整的n×n中间矩阵，而本方案仅计算n(n+1)/2个元素，浮点运算量减少约一半，大矩阵场景下性能提升明显。
• 避免高代价操作：无需对H进行平方根分解，省去了矩阵分解的额外开销。

内容的提问来源于stack exchange，提问作者SlamJam