In my project i work with covariance matrices and during many repetitive operation i have noticed that they gradually loose their symmetry. I put it down to numerical problems and when loss of symmetry exceed some level i try to restore matrix to proper form. I found in stack a simple method:
matrix = (matrix + matrix.transpose()) / 2.0
But when implemented i got strange result. Bottom triangle have proper values, but top one is different. Consider example
<code>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 << "nn";
matrix = (matrix + matrix.transpose()) / 2.0;
std::cout << std::fixed << matrix << std::scientific << "n"; // wrong upper triangle here
</code>
<code>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 << "nn";
matrix = (matrix + matrix.transpose()) / 2.0;
std::cout << std::fixed << matrix << std::scientific << "n"; // wrong upper triangle here
</code>
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 << "nn";
matrix = (matrix + matrix.transpose()) / 2.0;
std::cout << std::fixed << matrix << std::scientific << "n"; // wrong upper triangle here
When separate helper matrix used all if fine
<code> Eigen::Matrix< double, 9, 9 > helper = (matrix + matrix.transpose()) / 2.0;
matrix = helper;
std::cout << std::fixed << matrix << std::scientific << "n";
</code>
<code> Eigen::Matrix< double, 9, 9 > helper = (matrix + matrix.transpose()) / 2.0;
matrix = helper;
std::cout << std::fixed << matrix << std::scientific << "n";
</code>
Eigen::Matrix< double, 9, 9 > helper = (matrix + matrix.transpose()) / 2.0;
matrix = helper;
std::cout << std::fixed << matrix << std::scientific << "n";
Where is a problem?