Eigen 3.3 SVD.solve가 잘못된 값을 반환 함

Question

Eigen 3.2.7에서 3.3.4로 업데이트 한 후 JacobiSVD.solve에서 매우 잘못된 결과를 반환하는 문제가 발생했습니다. BDCSVD도 동일한 결과를 산출합니다.

문제는 다음 코드를 재현 할 수 있습니다 : 아이겐 3.3.4 계산 B4의

#include <Eigen/Eigen> 
#include <Eigen/SVD> 

int main(int argc, const char * argv[]) { 

    // M is calculated beforehand and does not change with different Eigen versions 
    Eigen::Matrix<double, 12, 12> M = Eigen::Matrix<double, 12, 12>::Zero(); 
    M(0,0) = 27; M(0,2) = 4.3625039999999995; M(0,3) = -9; M(0,5) = -2.1812519999999997; M(0,6) = -9; 
    M(1,1) = 27; M(1,2) = 3.2718720000000001; M(1,4) = -9; M(1,7) = -9; M(1,8) = -1.6359360000000001; 
    M(2,0) = 4.3625039999999995; M(2,1) = 3.2718720000000001; M(2,2) = 2.4780489612000003; 
    M(2,3) = -2.1812519999999997; M(2,5) = -0.82601632039999995; M(2,7) = -1.6359360000000001; 
    M(2,8) = -0.82601632039999995; M(3,0) = -9; M(3,2) = -2.1812519999999997; M(3,3) = 9; 
    M(4,1) = -9; M(4,4) = 9; M(5,0) = -2.1812519999999997; M(5,2) = -0.82601632039999995; 
    M(5,5) = 0.82601632039999995; M(6,0) = -9; M(6,6) = 9; M(7,1) = -9; M(7,2) = -1.6359360000000001; 
    M(7,7) = 9; M(8,1) = -1.6359360000000001; M(8,2) = -0.82601632039999995; M(8,8) = 0.82601632039999995; 

    Eigen::JacobiSVD<Eigen::MatrixXd> svd(M, Eigen::ComputeFullU); 
    Eigen::Matrix<double, 12, 12> ut = svd.matrixU(); 

    Eigen::Matrix<double, 6, 10> l_6x10; 
    Eigen::Matrix<double, 4, 6> dv0; 
    Eigen::Matrix<double, 4, 6> dv1; 
    Eigen::Matrix<double, 4, 6> dv2; 

    for(int i = 0; i < 4; i++) { 
     int a = 0, b = 1; 

     for(int j = 0; j < 6; j++) { 
      dv0(i, j) = ut(3 * a + 0, 11 - i) - ut(3 * b + 0, 11 - i); 
      dv1(i, j) = ut(3 * a + 1, 11 - i) - ut(3 * b + 1, 11 - i); 
      dv2(i, j) = ut(3 * a + 2, 11 - i) - ut(3 * b + 2, 11 - i); 

      b++; 
      if (b > 3) { 
       a++; 
       b = a + 1; 
      } 
     } 
    } 

    for(int i = 0; i < 6; i++) { 
     l_6x10(i,0) =  (dv0(0, i) * dv0(0, i) + dv1(0, i) * dv1(0, i) + dv2(0, i) * dv2(0, i)); 
     l_6x10(i,1) = 2.0f * (dv0(0, i) * dv0(1, i) + dv1(0, i) * dv1(1, i) + dv2(0, i) * dv2(1, i)); 
     l_6x10(i,2) =  (dv0(1, i) * dv0(1, i) + dv1(1, i) * dv1(1, i) + dv2(1, i) * dv2(1, i)); 
     l_6x10(i,3) = 2.0f * (dv0(0, i) * dv0(2, i) + dv1(0, i) * dv1(2, i) + dv2(0, i) * dv2(2, i)); 
     l_6x10(i,4) = 2.0f * (dv0(1, i) * dv0(2, i) + dv1(1, i) * dv1(2, i) + dv2(1, i) * dv2(2, i)); 
     l_6x10(i,5) =  (dv0(2, i) * dv0(2, i) + dv1(2, i) * dv1(2, i) + dv2(2, i) * dv2(2, i)); 
     l_6x10(i,6) = 2.0f * (dv0(0, i) * dv0(3, i) + dv1(0, i) * dv1(3, i) + dv2(0, i) * dv2(3, i)); 
     l_6x10(i,7) = 2.0f * (dv0(1, i) * dv0(3, i) + dv1(1, i) * dv1(3, i) + dv2(1, i) * dv2(3, i)); 
     l_6x10(i,8) = 2.0f * (dv0(2, i) * dv0(3, i) + dv1(2, i) * dv1(3, i) + dv2(2, i) * dv2(3, i)); 
     l_6x10(i,9) =  (dv0(3, i) * dv0(3, i) + dv1(3, i) * dv1(3, i) + dv2(3, i) * dv2(3, i)); 
    } 

    Eigen::Matrix<double, 6, 4> L_6x4; 

    for(int i = 0; i < 6; i++) { 
     L_6x4(i, 0) = l_6x10(i, 0); 
     L_6x4(i, 1) = l_6x10(i, 1); 
     L_6x4(i, 2) = l_6x10(i, 3); 
     L_6x4(i, 3) = l_6x10(i, 6); 
    } 

    // L_6x4 has the following values on 3.2.7 (everything else is 0): 
    // 
    // L_6x4(0,2) = 1; 
    // L_6x4(1,0) = 1; 
    // L_6x4(1,1) = 1; 
    // L_6x4(5,0) = -1.137432760287006; 
    // L_6x4(5,2) = -1.1374327602870071; 
    // L_6x4(5,3) = -1.1374327602870049; 
    // 
    // 
    // on 3.3.4 it has the following slightly different values: 
    // 
    // L_6x4(0,2) = 1; 
    // L_6x4(1,0) = 1; 
    // L_6x4(1,1) = 1; 
    // L_6x4(5,0) = -1.1374327602869998; 
    // L_6x4(5,2) = -1.1374327602870271; 
    // L_6x4(5,3) = -1.1374327602869889; 

    // Rho is calculated beforehand and does not change with different Eigen versions 
    Eigen::Matrix<double, 6, 1> Rho; 
    Rho << 0.25, 0.140625, 0, 0.390625, 0.25, 0.140625; 

    Eigen::JacobiSVD<Eigen::MatrixXd> l_6x4(L_6x4, Eigen::ComputeFullU | Eigen::ComputeFullV); 
    Eigen::Vector4d B4 = l_6x4.solve(Rho); 

    // The slight difference in L_6x4 apparently causes an insane difference in the result of solve. 
    // 
    // Eigen 3.2.7: B4={0.056766494562302421, 0, 0, -0.064568070601816963} 
    // Eigen 3.3.4: B4={-4773392957911.6992, 0, 0, -4196637484493.9165} 
    return 0; 
} 

값이 날이 몇 가지 이상한 메모리 정렬 문제의 SDK에서 버그가 수 있습니다 중 하나라고 생각하기 또는 잘하면이 프로그램의 버그.

-DEIGEN_DONT_ALIGN_STATICALLY 또는 -DEIGEN_DONT_VECTORIZE 및 -DEIGEN_DISABLE_UNALIGNED_ARRAY_ASSERT를 사용하여 컴파일하려고했습니다. 나는 또한이 매트릭스에서 Eigen :: DontAlign을 시도했다. 이들 중 어느 것도 결과에 아무런 영향을 미치지 않습니다.

는 NDK-r14b 그 소리의 C와 안드로이드 (4.4 및 5.1)에서 테스트 ++ _ 정적 및 Mac (맥 OS 시에라 10.12.6) : "그 소리 ++ -std = C++ (14) -g -O0 -I/EIGEN_PATH main.cpp -o main "

Answer 1

두 답변은 실제로 매우 유사하며 유사한 잔여를 산출합니다. 이는 행렬에 두 개의 0 특이 값이 있고 세 번째 값이 기계 정밀도에 매우 근접하기 때문입니다.

{0.056766494562302421, 0, 0, -0.064568070601816963} 

와 0.91의 잔류 : 따라서 반올림 오류에 따라,이 경우 당신은 솔루션의 나머지 3D 부분 공간 내에서 최소한의 표준 솔루션을 얻을 중 하나를 0으로 간주됩니다. 그렇지 않으면 의미로 간주됩니다, 그리고 당신은 (두 제로 특이 값에 해당하는) 작은 차원 부분 공간 내에서 최소한의 표준 솔루션을 얻을 : 0.89의 작은 잔류와

{-4773392957911.6992, 0, 0, -4196637484493.9165} 

. 따라서이 두 번째 해법은 좀 더 정확 합니다만, 더 작은 표준이지만 더 높은 잔차를 선호하는 경우 세 번째 특이 값을 0으로 간주하도록 임계 값을 조정할 수 있습니다. 이는 setThreshold을 통해 이루어집니다.

Eigen::JacobiSVD<Eigen::MatrixXd> l_6x4; 
l_6x4.setThreshold(1e-14); 
l_6x4.compute(L_6x4, Eigen::ComputeFullU | Eigen::ComputeFullV); 
Eigen::Vector4d B4 = l_6x4.solve(Rho);