Very large difference in results depending on the order of operations

Viewed 55

Depending on the sequence of operations of multiplication, the result for floating-point values is very different (for double = 0).

For demonstration, I use my matrix values (The result is clearer), but it will work for random initialization.

I understand that this is most likely due to float, but the difference is too big. this is especially noticeable for the first two elements of the final vector. It would not bother me if the difference between the elements of different vectors was 6-9 characters after the point, but the resulting difference greatly worsens the calculations.

    Eigen::MatrixXf m1 = Eigen::MatrixXf::Zero(12, 12);

    {
        m1(0, 0) = 4.5;
        m1(1, 1) = -4.5;
        m1(2, 3) = 4.5;
        m1(3, 4) = -4.5;
        m1(4, 6) = 4.5;
        m1(5, 7) = -4.5;
        m1(6, 0) = 4.5;
        m1(7, 1) = -4.5;
        m1(8, 9) = 4.5;
        m1(9, 10) = -4.5;
        m1(10, 6) = 4.5;
        m1(11, 7) = -4.5;
    }

    Eigen::MatrixXf m2 = Eigen::MatrixXf::Zero(18, 12);

    {
        m2(0, 0) = 119.87721252441406;
        m2(0, 1) = 119.8772201538086;
        m2(1, 0) = 119.87721252441406;
        m2(1, 1) = 119.8772201538086;
        m2(2, 0) = 119.87721252441406;
        m2(2, 1) = 119.8772201538086;
        m2(3, 2) = 119.87721252441406;
        m2(3, 3) = -119.8772201538086;
        m2(4, 2) = 119.87721252441406;
        m2(4, 3) = -119.8772201538086;
        m2(5, 2) = 119.87721252441406;
        m2(5, 3) = -119.8772201538086;
        m2(6, 4) = -119.87721252441406;
        m2(6, 5) = -119.8772201538086;
        m2(7, 4) = -119.87721252441406;
        m2(7, 5) = -119.8772201538086;
        m2(8, 4) = -119.87721252441406;
        m2(8, 5) = -119.8772201538086;
        m2(9, 6) = 119.87721252441406;
        m2(9, 7) = 119.8772201538086;
        m2(10, 6) = 119.87721252441406;
        m2(10, 7) = 119.8772201538086;
        m2(11, 6) = 119.87721252441406;
        m2(11, 7) = 119.8772201538086;
        m2(12, 8) = -119.87721252441406;
        m2(12, 9) = 119.8772201538086;
        m2(13, 8) = -119.87721252441406;
        m2(13, 9) = 119.8772201538086;
        m2(14, 8) = -119.87721252441406;
        m2(14, 9) = 119.8772201538086;
        m2(15, 10) = -119.87721252441406;
        m2(15, 11) = -119.8772201538086;
        m2(16, 10) = -119.87721252441406;
        m2(16, 11) = -119.8772201538086;
    }

    Eigen::MatrixXf m3 = Eigen::MatrixXf::Zero(18, 12);

    {
        m3(0, 0) = 535.1300048828125;
        m3(0, 1) = -535.1300659179688;
        m3(1, 0) = 535.1300048828125;
        m3(1, 1) = -535.1300659179688;
        m3(2, 0) = 535.1300048828125;
        m3(2, 1) = -535.1300659179688;

        m3(3, 3) = 535.1300048828125;
        m3(3, 4) = 535.1300659179688;
        m3(4, 3) = 535.1300048828125;
        m3(4, 4) = 535.1300659179688;
        m3(5, 3) = 535.1300048828125;
        m3(5, 4) = 535.1300659179688;

        m3(6, 6) = -535.1300048828125;
        m3(6, 7) = 535.1300659179688;
        m3(7, 6) = -535.1300048828125;
        m3(7, 7) = 535.1300659179688;
        m3(8, 6) = -535.1300048828125;
        m3(8, 7) = 535.1300659179688;

        m3(9, 0) = 535.1300048828125;
        m3(9, 1) = -535.1300659179688;
        m3(10, 0) = 535.1300048828125;
        m3(10, 1) = -535.1300659179688;
        m3(11, 0) = 535.1300048828125;
        m3(11, 1) = -535.1300659179688;

        m3(12, 9) = -535.1300048828125;
        m3(12, 10) = -535.1300659179688;
        m3(13, 9) = -535.1300048828125;
        m3(13, 10) = -535.1300659179688;
        m3(14, 9) = -535.1300048828125;
        m3(14, 10) = -535.1300659179688;

        m3(15, 6)= -535.1300048828125;
        m3(15, 7)= 535.1300659179688;
        m3(16, 6)= -535.1300048828125;
        m3(16, 7)= 535.1300659179688;
        m3(17, 6)= -535.1300048828125;
        m3(17, 7)= -535.1300048828125;
    }

    Eigen::MatrixXf ans1 = m1.transpose() * m2.transpose() * m3 * Eigen::VectorXf::Ones(12);

    Eigen::MatrixXf tmp1 = m3 * Eigen::VectorXf::Ones(12);
    Eigen::MatrixXf tmp2 = m2.transpose() * tmp1;
    Eigen::MatrixXf ans2 = m1.transpose() * tmp2;

    double error = (ans1 - ans2).norm();
error:
> 0.4130503535270691 
ans1:
> [0]  -0.125    
> [1]  0.375 
> [2]  0.0   
> [3]  1732047.25
> [4]  1732047.5 
> [5]  0.0   
> [6]  -0.125    
> [7]  0.25  
> [8]  0.0   
> [9]  1732047.25
> [10] 1732047.5 
> [11] 0.0 
ans2:
> [0]  -0.19755156338214874   
> [1]  0.19755156338214874
> [2]  0.0   
> [3]  1732047.25
> [4]  1732047.25 
> [5]  0.0   
> [6]  -0.16462630033493042    
> [7]  0.16462630033493042 
> [8]  0.0   
> [9]  1732047.25
> [10] 1732047.25 
> [11] 0.0 
0 Answers
Related