수치적으로 보다 정밀한 이차방정식의 해

영상처리에서 $2\times2$ 행렬의 고유값과 고유벡터를 수치적으로 구해야 하는 상황이 많이 생긴다. 이 경우 2차인 특성방정식을 풀어서 고유값을 구하게 되는데, 두 고유값의 차이가 매우 큰 경우 단순한 근의 공식을 적용한 결과는 수치적으로 부정확한 값을 줄 수 있다. 이를 피하는 방법을 간단히 살펴보자.

 

이차방정식 

$$x^2 - b x +c =0$$

의 두 해는 

\begin{align}x_1 = \frac{b + d}{2},\quad x_2 = \frac{b - d}{2},\qquad d= \sqrt{b^2 - 4c}\end{align}

로 표현된다.  고유값의 특성방정식으로 보면 계수 $b$는 두 고유값의 합, 판별식 $d$는 두 고유값의 차이다. 그런데 $b ^2 \gg c $ 판별식 $d$는 $b$와 매우 근접한 수치이므로 거의 비슷한 두 수의 뺄셈으로 표현된 $x_2$의 수치계산 결과는 매우 부정확해질 수 있다. 예로 

$$ x^2 - 10^8 x +1=0$$

의 해를 프로그램으로 구하면 

	double b = 1.0e8, c = 1.;
	double d = sqrt(b * b - 4. * c);
	double x1 = (b + d)/2., x2 = (b - d)/2.;
	TRACE("x1=%e, x2=%e", x1, x2);

프로그램을 돌리면 얻게 되는 결과는

$$ \tt  x1=1.000000e+008, \qquad x2=7.450581e-009$$

그런데 테일러 전개를 사용하면 $x_2$는 

$$x_2 = \frac{b}{2}\left( 1 - \sqrt{1-\frac{4c}{b^2}} \right)  \approx \frac{b}{2}\times \frac{2c}{b^2} = \frac{c}{b}=10^{-8}$$ 

이어야 한다. 따라서 $x_2$를 바로 근의 공식으로 구할 경우 25%정도의 오차가 생긴다. 이런 경우에는 근의 공식에 의한 직접적인 계산보다는 이차방정식의 특성을 이용하는 것이 수치적으로 더 정교한 답을 얻을 수 있다. 두 근의 곱이 $x_1 x_2 =c$이므로 $x_2$는 

$$x_2 = c / x_1\approx 10^{-8}$$

을 통해서 구하는 것이 더 수치적으로 정교하다.  일반적으로 두 근의 합  $b$가 양수인 경우는  $x_2$를, 음수인 경우는 $x_1$을 두 근의 곱을 이용해서 계산하는 것이 더 정교한 수치값을 얻을 수 있다.

// Find two roots of a quadratic eq: x * x - b * x + c = 0;
void Quadratic2(double b, double c, double& x1, double& x2) {
    double d = sqrt(b * b - 4. * c);
    x1 = (b + d)/ 2., x2 = (b - d)/ 2.;
    TRACE("x1=%e, x2=%e\r\n", x1, x2);
    if (b > 0) {
        if (x1 != 0) x2 = c / x1;
        TRACE("Improved: x1=%e, x2=%e\r\n", x1, x2);
    } else {
        if (x2 != 0) x1 = c / x2;
        TRACE("Improved: x1=%e, x2=%e\r\n", x1, x2);
    }
}

이 경우 결과는 훨씬 더 정교한 값을 보여준다.

$$\tt Improved: x1=1.000000e+008, \quad x2=1.000000e-008$$