2017-09-25 13 views
1

저는 지금 매트릭스 라이브러리를 사용하여 다양한 것을 계산하는 C++ 코드베이스를 만들고 있습니다. 그 중 하나는 행렬의 역수를 계산하는 것입니다. 그것을 달성하기 위해 가우스 제거를 사용합니다. 그러나 그 결과는 매우 부정확합니다. 순열 행렬을 원래 행렬에 곱하는 것은 항등식 행렬을 닫지 않습니다. 반대 인 경우gauss 제거를 사용하여 행렬의 역함수를 계산할 때 매우 높은 부정확함

/// \brief Take the inverse of the matrix. 
/// \return A new matrix which is the inverse of the current one. 
matrix<T, M, M> inverse() const 
{ 
    static_assert(M == N, "Inverse matrix is only defined for square matrices."); 

    // augmented the current matrix with the identiy matrix. 
    auto augmented = this->augment(matrix<T, M, M>::get_identity()); 

    for (std::size_t i = 0; i < M; i++) 
    { 
     // divide the current row by the diagonal element. 
     auto divisor = augmented[i][i]; 
     for (std::size_t j = 0; j < 2 * M; j++) 
     { 
      augmented[i][j] /= divisor; 
     } 

     // For each element in the column of the diagonal element that is currently selected 
     // set all element in that column to 0 except the diagonal element by using the currently selected row diagonal element. 

     for (std::size_t j = 0; j < M; j++) 
     { 
      if (i == j) 
      { 
       continue; 
      } 

      auto multiplier = augmented[j][i]; 
      for (std::size_t k = 0; k < 2 * M; k++) 
      { 
       augmented[j][k] -= multiplier * augmented[i][k]; 
      } 
     } 
    } 

    // Slice of the the new identity matrix on the left side. 
    return augmented.template slice<0, M, M, M>(); 
} 

지금은 단위 테스트 시험을 한 다음은

의 역을 계산하는 데 사용되는 코드는, 매트릭스가 수치 형에 주형하고, 행과 열이 사전 계산 된 값을 사용하여 정정하십시오. 나는 하나의 3x3과 하나의 4x4의 두 행렬을 시도한다. 이 웹 사이트를 사용하여 역수를 계산했습니다 : https://matrix.reshish.com/ 그리고 어느 정도 일치합니다. 단위 테스트가 성공했기 때문입니다. 그러나 일단 원래의 행렬을 계산하면 역행렬은 항등 행렬과 닮았습니다. 아래 코드의 주석을 참조하십시오.

BOOST_AUTO_TEST_CASE(matrix_inverse) 
{ 
    auto m1 = matrix<double, 3, 3>({ 
     {7, 8, 9}, 
     {10, 11, 12}, 
     {13, 14, 15} 
    }); 

    auto inverse_result1 = matrix<double,3, 3>({ 
     {264917625139441.28, -529835250278885.3, 264917625139443.47}, 
     {-529835250278883.75, 1059670500557768, -529835250278884.1}, 
     {264917625139442.4, -529835250278882.94, 264917625139440.94} 
    }); 

    auto m2 = matrix<double, 4, 4>({ 
     {7, 8, 9, 23}, 
     {10, 11, 12, 81}, 
     {13, 14, 15, 11}, 
     {1, 73, 42, 65} 
    }); 

    auto inverse_result2 = matrix<double, 4, 4>({ 
     {-0.928094660194201, 0.21541262135922956, 0.4117111650485529, -0.009708737864078209}, 
     {-0.9641231796116679, 0.20979975728155775, 0.3562651699029188, 0.019417475728154842}, 
     {1.7099261731391882, -0.39396237864078376, -0.6169346682848 , -0.009708737864076772 }, 
     {-0.007812499999999244, 0.01562499999999983, -0.007812500000000278, 0} 
    }); 


    // std::cout << (m1.inverse() * m1) << std::endl; 
    // results in 
    // 0.500000000  1.000000000  -0.500000000 
    // 1.000000000  0.000000000  0.500000000 
    // 0.500000000  -1.000000000 1.000000000 
    // std::cout << (m2.inverse() * m2) << std::endl; 
    // results in 
    // 0.396541262  -0.646237864 -0.689016990 -2.162317961 
    // 1.206917476  2.292475728  1.378033981  3.324635922 
    // -0.884708738 -0.958737864 -0.032766990 -3.756067961 
    // -0.000000000 -0.000000000 -0.000000000 1.000000000 

    BOOST_REQUIRE_MESSAGE(
     m1.inverse().fuzzy_equal(inverse_result1, 0.1) == true, 
     "3x3 inverse is not the expected result." 
    ); 

    BOOST_REQUIRE_MESSAGE(
     m2.inverse().fuzzy_equal(inverse_result2, 0.1) == true, 
     "4x4 inverse is not the expected result." 
    ); 
} 

나는 재치가있다. 나는 수학에서 수학을 전문으로하는 사람은 아니지만, 직업에서이 모든 것을 배워야했기 때문에 이것은 정말로 나를 괴롭 히고있다. 역함수의 위치 https://codeshare.io/johnsmith

선로 (404)는 :

전체 코드 행렬 클래스에서 확인할 수있다.

도움을 주시면 감사하겠습니다.

+2

일부 매트릭스에는 역 매트릭스가 없습니다. 처음 3x3 예제와 같습니다. 그들을 위해, 프로그램은 말도 안되는 결과를 만들어 낼 것입니다. – anatolyg

+2

행렬의 행렬식이 0이 아닌지 확인하십시오 (즉, 행렬이 비 고유이며 반전 될 수 있음). 0이 아닌 0에 매우 가까운 행렬식 인 경우에도 조심해야합니다 (조건 번호에 대해 읽음). https://math.stackexchange.com/questions/1622610/when-is-inverting-a-matrix-numerically-unstable – vsoftco

+0

@vsoftco 두 경우 모두 행렬식은 0이 아닙니다. –

답변

5

주석에 이미 설정된대로 관심있는 매트릭스는 단수이므로 역수는 없습니다.

테스트 결과 코드에서 이미 첫 번째 문제가 발견되었습니다.이 경우 제대로 처리되지 않으며 오류가 발생하지 않습니다.

큰 문제는 감지하기 쉽지 않다는 것입니다. 반올림 오류로 인한 오류가없는 경우 조각의 케이크가됩니다. 제수가 0이 아닌지 테스트하십시오. 그러나 부동 연산에 반올림 오류가 있으므로 제수가 0이 아닌 매우 작은 숫자가됩니다.

반올림 오류로 인해이 0이 아닌 값이 있는지 또는 행렬이 단수 (단수가 아님) 근처에 있는지 알 수있는 방법이 없습니다. 그러나 매트릭스가 단수에 가깝다면 상태가 좋지 않아 결과를 신뢰할 수 없습니다.

이상적으로 알고리즘은 역행렬을 계산할뿐만 아니라 원래 행렬의 조건도 추정해야하므로 호출자는 나쁜 조건에 대해 반응 할 수 있습니다.

아마도 이런 종류의 계산에는 잘 알려지고 잘 테스트 된 라이브러리를 사용하는 것이 좋습니다. 고려해야 할 사항이 많고 잘못 수행 할 수있는 사항이 있습니다.