Reply to topic

Numerical instability in .inverse()?

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
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
User avatar ggael
Moderator
Posts
2193
Karma
15
OS
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
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!

 
Reply to topic

Bookmarks



Who is online

Registered users: AGB, Azzart, Baidu [Spider], Bing [Bot], Bingy, colomar, Exabot [Bot], garthecho, Google [Bot], google01103, GreatEmerald, jensreuterberg, jstaniek, kainz.a, koriun, Kver, La Ninje, lazyit, Majestic-12 [Bot], MSN [Bot], Nuc!eoN, psbot [Picsearch], rulet111, scummos, SecretCode, Tuukka, undoshmulhe, Uri_Herrera, Yahoo [Bot]