Numerical instability in .inverse()?

Board index

Page 1 of 1 (4 posts)

Tags:

damien_d Registered Member Posts 15 Karma 0	Numerical instability in .inverse()? Wed Jun 12, 2013 12:49 am Dear all, I am trying to get my head around a result that I have seen in my code. In short, I have a snippet of code: Code: Select all template <typename Derived1, typename Derived2, typename Derived3, typename Derived4, typename Derived5> void Estimators::KalmanFilter::PartiallyObservedObservationModel( const Eigen::MatrixBase<Derived1>& mu_obs_p, const Eigen::MatrixBase<Derived2>& mu_un_p, const Eigen::MatrixBase<Derived3>& C_obs_obs_p, // Top left of covariance matrix const Eigen::MatrixBase<Derived4>& C_un_obs_p, // Bottom Left of covariance matrix const Eigen::MatrixBase<Derived5>& C_un_un_p, // Bottom right of covariance matrix const Eigen::MatrixBase<Derived1>& mu_obs_e, Eigen::MatrixBase<Derived2>& mu_un_e, const Eigen::MatrixBase<Derived3>& C_obs_obs_e, // Top left of covariance matrix Eigen::MatrixBase<Derived4>& C_un_obs_e, // Bottom Left of covariance matrix Eigen::MatrixBase<Derived5>& C_un_un_e) // Bottom right of covariance matrix { // [1] Un-numbered equation after (5). const Eigen::MatrixBase<Derived4> L = C_un_obs_p * C_obs_obs_p.inverse(); std::cout << "C_un_obs_p = " << std::endl << C_un_obs_p << std::endl << std::endl; std::cout << "C_obs_obs_p = " << std::endl << C_obs_obs_p << std::endl << std::endl; std::cout << "L = " << std::endl << L << std::endl << std::endl; Which produces the following result: Code: Select all C_un_obs_p = -0.000160881 0.00469207 0.000465655 1.25005 2.27872e-05 -2.43644e-05 0 0 0 -0.00611925 0.00348598 0.000684872 1.36718e-05 1.25008 1.16715e-05 0 0 0 -0.0045807 -0.00489111 -0.00021718 -2.43644e-05 1.16715e-05 1.25007 0 0 0 0 0 0 1.71495e-06 1.44875e-06 -1.99801e-06 0 0 0 0 0 0 3.02184e-07 2.29146e-06 1.9209e-06 0 0 0 0 0 0 -2.44941e-06 1.29704e-06 -1.16192e-06 0 0 0 -1.52309e-07 7.16848e-11 2.35168e-11 2.01558e-11 7.53078e-10 5.63061e-10 0 0 0 -7.16916e-11 -1.52309e-07 -3.27793e-11 -5.77202e-10 -4.2831e-10 6.02072e-10 0 0 0 -2.34962e-11 3.27941e-11 -1.52309e-07 -5.70538e-11 -8.39552e-11 2.66576e-11 0 0 0 C_obs_obs_p = 0.0171348 -2.09192e-12 -1.44083e-11 -4.02202e-06 -0.000152981 -0.000114518 0 0 0 -2.09192e-12 0.0171348 2.25817e-11 0.000117302 8.71494e-05 -0.000122278 0 0 0 -1.44083e-11 2.25817e-11 0.0171348 1.16414e-05 1.71218e-05 -5.42951e-06 0 0 0 -4.02202e-06 0.000117302 1.16414e-05 25.0625 4.09704e-07 -5.56414e-07 0 0 0 -0.000152981 8.71494e-05 1.71218e-05 4.09704e-07 25.0625 2.64503e-07 0 0 0 -0.000114518 -0.000122278 -5.42951e-06 -5.56414e-07 2.64503e-07 25.0625 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 L = 1.21892e-314 2.84789e-306 -7.48068e+292 3.0138e-322 2.84789e-306 2.84789e-306 3.11261e-322 0 0 1.96759e-314 2.84789e-306 2.84789e-306 -5.25495e-79 0 2.84789e-306 -8870.89 0 0 2.84789e-306 2.84789e-306 -4.31733e-297 4.52231e-66 1.96724e-314 2.27665e-320 -3.86837e+191 3.32632e-212 0 2.84789e-306 0 9.33784e-322 -6.07507e-35 0 2.27665e-320 4.80136e+64 9.48606e-322 0 2.84789e-306 4.41086e+92 -3.75209e-224 4.41086e+92 0 2.27665e-320 9.43665e-322 0 0 1.96726e-314 2.56367e+18 3.11261e-322 -3.75209e-224 2.56367e+18 -8870.89 0 0 0 2.84789e-306 1.96831e-314 -3.75209e-224 2.84385e+31 2.84385e+31 -6.61981e-70 -8870.89 0 0 6.14962e-315 4.41086e+92 3.11261e-322 2.84789e-306 2.84789e-306 -1.09718e+146 2.17005e-196 0 0 2.84789e-306 -1.47242e+257 -6.07507e-35 -5.07662e+261 2.84789e-306 3.88938e-308 9.43665e-322 0 0 Note the extremely large values I have highlighted. When using the same in MATLAB: Code: Select all >> C_un_obs_p * inv(C_obs_obs_p) ans = -0.0094 0.2735 0.0271 0.0499 -0.0000 0.0000 0 0 0 -0.3567 0.2032 0.0399 -0.0000 0.0499 -0.0000 0 0 0 -0.2670 -0.2851 -0.0127 0.0000 -0.0000 0.0499 0 0 0 -0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 0 0 0 0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 0 0 0 0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 0 0 0 -0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 0 0 0 -0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 0 0 0 -0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 0 0 0 >> C_un_obs_p / C_obs_obs_p ans = -0.0094 0.2735 0.0271 0.0499 -0.0000 0.0000 0 0 0 -0.3567 0.2032 0.0399 -0.0000 0.0499 -0.0000 0 0 0 -0.2670 -0.2851 -0.0127 0.0000 -0.0000 0.0499 0 0 0 -0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 0 0 0 0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 0 0 0 0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 0 0 0 -0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 0 0 0 -0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 0 0 0 -0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 0 0 0 Which is about what I expect. Both are passed into the function as fixed-size matricies: Code: Select all Estimators::KalmanFilter::PartiallyObservedObservationModel( mu_do_kminus, // Directly observed states prior to update mu_io_kminus, // Indirectly observed states prior to update Pdd_kminus, // Top left of covariance matrix Pid_kminus, // Bottom Left of covariance matrix Pii_kminus, // Bottom right of covariance matrix mu_do_kplus, mu_io_kplus, Pdd_kplus, // Top left of covariance matrix Pii_kplus, // Bottom Left of covariance matrix Pid_kplus); // Bottom right of covariance matrix Where: Code: Select all `Eigen::Matrix<double, 9, 9> Pdd_kminus; Eigen::Matrix<double, 9, 9> Pid_kminus; Eigen::Matrix<double, 9, 9> Pii_kminus;` So I'm not sure why these values are produced. Is anyone able to shed light as to what may be going wrong? -- Damien
damien_d Registered Member Posts 15 Karma 0	Re: Numerical instability in .inverse()? Wed Jun 12, 2013 1:08 am I am convinced it is not the inverse itself. I have just added some instrumentation with the following result: Code: Select all C_un_obs_p = -8.98667e-06 0.00469869 8.9643e-05 1.25005 2.4617e-05 -2.31177e-05 0 0 0 -0.00582331 0.00403296 5.81367e-06 1.55016e-05 1.25008 1.15399e-05 0 0 0 -0.00508335 -0.00462559 -0.000150529 -2.31177e-05 1.15399e-05 1.25007 0 0 0 0 0 0 1.85507e-06 1.55266e-06 -1.78323e-06 0 0 0 0 0 0 4.9899e-08 2.24034e-06 2.00258e-06 0 0 0 0 0 0 -2.36393e-06 1.26573e-06 -1.3571e-06 0 0 0 -1.52309e-07 4.34401e-11 -6.21741e-11 1.30037e-12 7.16505e-10 6.25072e-10 0 0 0 -4.34658e-11 -1.52309e-07 4.71197e-11 -5.78015e-10 -4.95767e-10 5.69238e-10 0 0 0 6.21561e-11 -4.71433e-11 -1.52309e-07 -1.12157e-11 -1.38164e-12 1.8449e-11 0 0 0 C_obs_obs_p = 0.0171348 -2.15168e-12 -1.93664e-11 -2.24659e-07 -0.000145583 -0.000127084 0 0 0 -2.15168e-12 0.0171348 4.27406e-13 0.000117467 0.000100824 -0.00011564 0 0 0 -1.93664e-11 4.27406e-13 0.0171348 2.24108e-06 1.45343e-07 -3.76323e-06 0 0 0 -2.24659e-07 0.000117467 2.24108e-06 25.0625 4.55449e-07 -5.25246e-07 0 0 0 -0.000145583 0.000100824 1.45343e-07 4.55449e-07 25.0625 2.61448e-07 0 0 0 -0.000127084 -0.00011564 -3.76323e-06 -5.25246e-07 2.61448e-07 25.0625 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 C_obs_obs_p.inv() = 58.3609 6.14356e-09 1.28011e-07 5.23144e-07 0.000339006 0.000295929 0 0 0 6.14356e-09 58.3609 9.54527e-08 -0.000273536 -0.00023478 0.00026928 0 0 0 1.28011e-07 9.54527e-08 58.3609 -5.21861e-06 -3.38448e-07 8.7631e-06 0 0 0 5.23144e-07 -0.000273536 -5.21861e-06 0.0399002 3.7839e-10 -4.24035e-10 0 0 0 0.000339006 -0.00023478 -3.38448e-07 3.7839e-10 0.0399002 2.19415e-10 0 0 0 0.000295929 0.00026928 8.7631e-06 -4.24035e-10 2.19415e-10 0.0399002 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 L = 1.27124e-314 2.84796e-306 2.57463e+18 2.84796e-306 3.06321e-322 0 0 0 0 5.48681e-314 2.84796e-306 2.84796e-306 2.84796e-306 -8870.89 0 0 0 0 2.84796e-306 2.84796e-306 4.40288e+92 2.27665e-320 -2.18508e-247 3.32632e-212 0 0 0 2.84796e-306 0 2.57464e+18 2.27665e-320 7.83061e+64 9.48606e-322 0 0 0 5.48648e-314 0 2.85695e+31 2.27665e-320 9.43665e-322 0 0 0 0 5.48646e-314 2.84796e-306 9.53547e-322 -8870.89 0 0 0 0 0 2.84796e-306 3.58518e-26 0 -1.74107e-67 -8870.89 0 0 0 0 6.16589e-315 2.84796e-306 2.84796e-306 -2.12445e+175 -8.23733e+240 0 0 0 0 2.84796e-306 4.40288e+92 2.84796e-306 -7.57956e-16 9.38725e-322 0 0 0 0 Which is "pretty close" to MATLAB's inv() Code: Select all >> inv(C_obs_obs_p) - C_obs_obs_p_inv ans = 1.0e-03 * -0.1320 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0 0 0 0.0000 -0.1320 -0.0000 0.0000 0.0000 0.0000 0 0 0 -0.0000 -0.0000 -0.1371 0.0000 0.0000 -0.0000 0 0 0 -0.0000 0.0000 0.0000 0.0001 -0.0000 0.0000 0 0 0 -0.0000 0.0000 0.0000 -0.0000 0.0001 0.0000 0 0 0 -0.0000 0.0000 -0.0000 0.0000 0.0000 0.0001 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 But nowhere near the actual result (note the factor of 10^240): Code: Select all >> (C_un_obs_p * C_obs_obs_p_inv) - L ans = 1.0e+240 * -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 0 0 0 -0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 0 0 0 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 0.0000 0 0 0 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0 0 0 0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 0 0 0 0.0000 -0.0000 0.0000 0.0000 0.0000 -0.0000 0 0 0 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 0 0 0 -0.0000 -0.0000 0.0000 0.0000 8.2373 -0.0000 0 0 0 0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 0 0 0 So, is there something wrong with the multiplcation and/or the assignment? -- Damien
ggael Moderator Posts 3447 Karma 19 OS	Re: Numerical instability in .inverse()? Wed Jun 12, 2013 8:49 am I'm surprised this program even compile or run without assertions. The statement: Code: Select all `Eigen::MatrixBase<Derived4> L = ...` is invalid. A MatrixBase<> has no storage, you cannot create such an object. Use: Code: Select all `typename Derived4::PlainObject L =` Moreover, you'll get faster and more accurate result without explicitly inverting the matrix, e.g.: Code: Select all `L.transpose() = C_obs_obs_p.transpose().lu().solve(C_un_obs_p.transpose());` In the future, Eigen will do this rewriting for you. Currently avoid .inverse() unless you really want it!
damien_d Registered Member Posts 15 Karma 0	Re: Numerical instability in .inverse()? Thu Jun 13, 2013 12:32 am ggael wrote:I'm surprised this program even compile or run without assertions. So was I ggael wrote:Moreover, you'll get faster and more accurate result without explicitly inverting the matrix, e.g.: Code: Select all `L.transpose() = C_obs_obs_p.transpose().lu().solve(C_un_obs_p.transpose());` In the future, Eigen will do this rewriting for you. Currently avoid .inverse() unless you really want it! Agreed, but I was trying to find the problem that has now been tracked to the storage... thanks!