Getting Double Precision(MatrixXd) answers using Single Precision(MatrixXf) on Eigen C++

Question

I have coded a model predictive controller on a pixhawk hardware to control the attitude (angles) of the quadcopter. I exchanged messages with one of the developers of pixhawk and he advised me to use single precision. My code is in double precision.

Before this, I tested a numerical problem with MatrixXd (double precision) using Eigen C++ Library and my code was able to get the same answer. I used the ldlt() Cholesky solver for a dense linear systems (all the other solver methods resulted in the wrong answer).

I replaced all MatrixXd with MatrixXf and double with float in order to solve the problem with single precision but I was unable to get the same answer. So I want to get some insight as to why I am not getting the same answer achieved by MatrixXd when I use MatrixXf.

Below is the relevant section. The y variable at the end produces the -nan(ind); -nan(ind) as it's solution when I declare my variables and matrices with MatrixXf but when I use MatrixXd I get the desired solution of 1; 0.9999.

// Hildreth's Quadratic Programming Loop

for (int i = 0; i < 3; i++)//i < r_cols - 1; i++)
{
    MatrixXf F = -2 * (H.transpose())*(Rs*r.col(i) - P*Xf);
    MatrixXf d = dd + dupast*uin;
    MatrixXf DeltaU = QPhild(E, F, CC, d);

    MatrixXf DeltaU_1(4, 2);
    DeltaU_1 << DeltaU(0, 0), DeltaU(1, 0),
        DeltaU(2, 0), DeltaU(3, 0),
        DeltaU(4, 0), DeltaU(5, 0),
        DeltaU(6, 0), DeltaU(7, 0);

    MatrixXf deltau = DeltaU_1.row(0);
    MatrixXf deltau_tran = deltau.transpose();
    u = u + deltau_tran;

    // Process
    x.col(i + 1) = Ad*x.col(i) + Bd*u;
    y = Cd*x.col(i + 1) + dist.col(i);

    // Model
    xh.col(i + 1) = A*xh.col(i) + B*deltau_tran + L*(y - C*xh.col(i));
    yh = C*xh.col(i + 1);

    Xf << x.col(i + 1) - x.col(i),
        y;
  }

 cout << y << endl << endl;

Below is the QPhild function:

MatrixXf QPhild(MatrixXf E, MatrixXf F, MatrixXf CC, MatrixXf d)
{
   MatrixXf CC_trans = CC.transpose();
   MatrixXf T = CC*(E.ldlt().solve(CC_trans));
   MatrixXf K = (CC*(E.ldlt().solve(F)) + d);

   int k_row = K.rows();
   int k_col = K.cols();

   MatrixXf lambda(k_row, k_col);
   lambda.setZero(k_row, k_col);

   MatrixXf al(0, 0);
   al.setConstant(10.0f);

   for (int km = 0; km < 40; km++)
   {
       MatrixXf lambda_p = lambda;

      // loop to determine lambda values for respective iterations
      for (int i = 0; i < k_row; i++)
      {
        MatrixXf t1 = T.row(i)*lambda;

        float t2 = T(i, i)*lambda(i, 0);
        float w = t1(0, 0) - t2;

        w = w + K(i, 0);
        float la = -w / T(i, i);

        if (la < 0.0f)  lambda(i, 0) = 0.0f;
        else lambda(i, 0) = la;
       }
       al = (lambda - lambda_p).transpose() * (lambda - lambda_p);

       float all = al(0, 0);
       float tol = 0.0000001f;

       if (all < tol) break;
   }

   MatrixXf DeltaU = -E.ldlt().solve(F) - (E.ldlt().solve(CC_trans))*lambda;

   return DeltaU;
}

I know that the main problem is from the function above because I put various cout lines to check outputs of various matrices within it. They start off at 0 then go to inf then to nan.

EDIT

Right now I am using the LLT decomposition and not LDLT (although they both give me the same answers). Regardless, I am posting the values of Matrix E in both double and float, with the corresponding singular values.

Matrix E in double:

1.84805e+12 1.65144e+12   7.557e+11 6.73531e+11 3.08645e+11 2.73966e+11 1.25821e+11 1.10981e+11
1.65144e+12 1.47576e+12 6.75306e+11 6.01881e+11 2.75811e+11 2.44823e+11 1.12436e+11 9.91757e+10
  7.557e+11 6.75306e+11  3.0902e+11  2.7542e+11 1.26211e+11  1.1203e+11 5.14507e+10 4.53824e+10
6.73531e+11 6.01881e+11  2.7542e+11 2.45475e+11 1.12488e+11 9.98504e+10 4.58567e+10 4.04488e+10
3.08645e+11 2.75811e+11 1.26211e+11 1.12488e+11 5.15474e+10  4.5756e+10 2.10137e+10 1.85354e+10
2.73966e+11 2.44823e+11  1.1203e+11 9.98504e+10  4.5756e+10 4.06159e+10 1.86529e+10 1.64534e+10
1.25821e+11 1.12436e+11 5.14507e+10 4.58567e+10 2.10137e+10 1.86529e+10 8.56643e+09  7.5562e+09
1.10981e+11 9.91757e+10 4.53824e+10 4.04488e+10 1.85354e+10 1.64534e+10  7.5562e+09 6.66534e+09

~~I cannot put all the numbers on a single line but the space between the rows is used to distinguish between different rows.~~ Singular values of E in double:

3.98569e+12 5.24887e+06  1363.09  174.56  166.311  159.098  58.9402  54.5173

Matrix E in float:

1.84805e+12 1.65144e+12   7.557e+11 6.73531e+11 3.08645e+11 2.73966e+11 1.25821e+11 1.10981e+11
1.65144e+12 1.47576e+12 6.75307e+11 6.01881e+11 2.75811e+11 2.44823e+11 1.12436e+11 9.91757e+10
  7.557e+11 6.75307e+11  3.0902e+11  2.7542e+11 1.26211e+11  1.1203e+11 5.14507e+10 4.53824e+10
6.73531e+11 6.01881e+11  2.7542e+11 2.45475e+11 1.12488e+11 9.98505e+10 4.58567e+10 4.04488e+10
3.08645e+11 2.75811e+11 1.26211e+11 1.12488e+11 5.15474e+10  4.5756e+10 2.10137e+10 1.85354e+10
2.73966e+11 2.44823e+11  1.1203e+11 9.98505e+10  4.5756e+10 4.06159e+10 1.86529e+10 1.64534e+10
1.25821e+11 1.12436e+11 5.14507e+10 4.58567e+10 2.10137e+10 1.86529e+10 8.56643e+09  7.5562e+09
1.10981e+11 9.91757e+10 4.53824e+10 4.04488e+10 1.85354e+10 1.64534e+10  7.5562e+09 6.66534e+09

And the float singular values are:

3.98569e+12 5.08741e+06  62753.4  58133.2  26340.8  20529.1 15839.4 1050.96

Without looking into any details of your code, the problem likely is that you need to increase `tol` when you work with `float` instead of `double`. — chtz, Aug 26 '17 at 17:41
I just did that, but sadly there was no change, I still got the same answer as before — Chinedu Amadi, Aug 26 '17 at 18:29
Please show where exactly you get NaNs and what the values of your matrices are (i.e., provide a [mcve]!) Overall, you do lots of sub-optimal stuff, like decomposing `E` multiple times, passing matrices by value. Is your question about sparse to dense conversion related, to the `double` vs `float` problem? Otherwise, ask a separate question for that (and reduce the question to the actual problem, don't spam too much irrelevant context) — chtz, Aug 26 '17 at 18:41
I have made some edits to the question, please let me know if it is clearer now. I have taken notice of your suggestions, I will pass my matrices by reference instead. This is my first time programming in linear algebra problems in C++, how do you propose I go about correcting multiple decompositions of `E`? — Chinedu Amadi, Aug 26 '17 at 19:25
Your code is still far away from being a MCVE, neither is it minimal, nor is it compilable (as it is). Reduce your code to the part which generates the `nan` -- I don't want to debug through your code for you. Regarding decomposing `E` just once: Declare a decomposition `const Eigen::LDLT E_ldlt(E);` once and use it's `.solve()` method wherever needed. — chtz, Aug 26 '17 at 19:38
Okay thank you. I have made it minimal but I cannot make it compilable as there are functions and variables that are needed and putting all that is essentially putting everything..about 500+ lines of code. I started trying the other methods again `llt()` gives me real numbers but they are far off about `2.6; 3.15` — Chinedu Amadi, Aug 26 '17 at 21:55
As commented elsewhere: Hard to help without knowing the contents of the matrix E. If only LDLT is giving you a satisfactory result when using double, then this means there is something fishy about it. Since it seems to be small you could past it's content here (with full accuracy, both using double and float), and also print its singular-values: `cout << E.jacobiSvd().singularValues().transpose() << endl;`, again both in double and single precision. — ggael, Aug 30 '17 at 11:39
Hello Ggael, I have made an edit to the post and added the values you requested — Chinedu Amadi, Aug 30 '17 at 20:36
Your matrix is near-singular when using single precision. You should therefore use a rank-revealing decomposition, like `FullPivLU`, `ColPivHouseholderQR`, or `SelfAdjointEigenSolver`. See here for an overview of available decompositions: http://eigen.tuxfamily.org/dox/group__TopicLinearAlgebraDecompositions.html — chtz, Sep 02 '17 at 09:19
Hello Chtz, I am not completely sure of what you wanted me to do, but this is what I did. I first of all declared this: `const ColPivHouseholderQR ColofE(E);` and used `ColofE` anytime I needed to use the `solve()` method. The answer i got fro this method is `2.65872; 3.14192` which is not close to the desired `1.00;1.00`. I tried doing the QPhild solution in `MatrixXd` then casting to `MatrixXf` but after a while it goes to infinity then `nan`. — Chinedu Amadi, Sep 03 '17 at 20:25

Getting Double Precision(MatrixXd) answers using Single Precision(MatrixXf) on Eigen C++

0 Answers0