Очень большая разница в результатах в зависимости от порядка операций - PullRequest
Новая
       94

Очень большая разница в результатах в зависимости от порядка операций

2 голосов
/ 13 июля 2020

В зависимости от последовательности операций умножения результат для значений с плавающей запятой сильно отличается (для double = 0).

Для демонстрации я использую свои матричные значения (результат более понятен), но он будет работать для случайной инициализации.

Я понимаю, что это, скорее всего, из-за float, но разница слишком велика. особенно это заметно для первых двух элементов конечного вектора. Меня бы не беспокоило, если бы разница между элементами разных векторов была на 6-9 знаков после точки, но полученная разница сильно ухудшает вычисления. 

    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
...