使用Eigen修正协方差矩阵对称性时出现异常结果的问题

问题场景 #

处理协方差矩阵时，因数值计算误差导致矩阵逐渐失去对称性，尝试用matrix = (matrix + matrix.transpose()) / 2.0原地恢复对称性时，出现下三角值正确、上三角值错误的异常；但先将结果存入辅助矩阵再赋值给原矩阵时，结果完全正常。

异常代码示例 #

Eigen::Matrix< double, 9, 9 > matrix;
matrix << 40559.778825, -31153.274570 , -2446.164237 , 17250.180423, -13036.183685 ,  -614.749476 ,  4145.292673,  -2929.325056 ,   -76.389702,
        -31163.652992 , 24642.034288 , -3155.954773 ,-13024.935948 , 10606.170622 ,  -793.926649 , -2912.511277  , 2669.226555 ,   -98.352932,
        -2446.293563 , -3155.849492,  81838.610505,   -618.601965 ,  -798.172922 , 20679.490319 ,   -80.854530 ,  -103.928341  , 2647.708959,
        17248.482322, -13019.753998 ,  -618.571123,  12761.733661 , -9040.721303,   -280.843475  , 5353.243186 , -3137.980275  ,  -66.415783,
        -13042.257916 , 10607.886484 ,  -798.195886 , -9042.370514  , 8146.254203 ,  -362.980644  ,-3133.173126 ,  3755.211086  ,  -85.807600,
        -614.817862 ,  -793.873126 , 20679.490338 ,  -280.857733  , -362.970017 ,  9637.963651 ,   -67.965019  ,  -87.712623 ,  2515.917137,
        4145.394398,  -2912.244650 ,   -80.851248  , 5353.380853,  -3133.239131  ,  -67.964682  , 4019.412600 , -1611.914593  ,  -29.487958,
        -2929.547642 ,  2669.130747 ,  -103.929154 , -3137.845462,   3755.074471 ,   -87.712468,  -1611.816127  , 3196.926413 ,   -38.111306,
        -76.393072 ,   -98.350422 ,  2647.708977,    -66.413746,    -85.809278 ,  2515.917141 ,   -29.487009  ,  -38.112091,   1382.242447;

std::cout << std::fixed << matrix << std::scientific << "\n\n";
matrix = (matrix + matrix.transpose()) / 2.0;
std::cout << std::fixed << matrix << std::scientific << "\n"; // 上三角值错误

正常代码示例 #

Eigen::Matrix< double, 9, 9 > matrix;
// 初始化代码同上述异常示例

std::cout << std::fixed << matrix << std::scientific << "\n\n";
Eigen::Matrix< double, 9, 9 > helper = (matrix + matrix.transpose()) / 2.0;
matrix = helper;
std::cout << std::fixed << matrix << std::scientific << "\n"; // 结果正常

原因分析 #

这是Eigen的原地赋值内存覆盖问题导致的：

• 当执行matrix = (matrix + matrix.transpose()) / 2.0时，Eigen会尝试进行原地优化，在计算过程中直接修改原矩阵的元素。
• 计算顺序上，Eigen可能先更新下三角区域的元素，后续计算上三角区域时，需要用到的原矩阵上三角元素已经被之前的修改覆盖（转置操作会引用原矩阵元素，但此时原矩阵部分元素已被更新），最终导致上三角区域计算错误。
• 使用辅助矩阵时，(matrix + matrix.transpose()) / 2.0的计算完全基于原矩阵的原始值，所有元素计算完成后再赋值给原矩阵，避免了内存覆盖问题，因此结果正确。

额外建议 #

如果想避免使用辅助矩阵，也可以调用Eigen专门提供的对称性恢复函数：

matrix.selfadjointView<Eigen::Lower>().rankUpdate(matrix, 0.5, 0.5);

或者明确要求Eigen先计算临时矩阵再赋值：

matrix = Eigen::Matrix<double,9,9>(matrix + matrix.transpose()) / 2.0;

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