Bairstow's method
안녕하세요?
오늘은 다항식의 근사해를 찾아내는 수치해석 알고리즘인 Bairstow’s method에 관해 알아보겠습니다.
다항식의 근사해 찾기
주어진 다항식의 해를 찾는 방법은 어떤 것들이 있을까요?
가장 먼저 생각해 볼 수 있는 것은 근의 공식입니다. 예를 들면 $ax^2 + bx + c = 0$이라는 다항식의 해가 $\frac{-b \pm \sqrt{b^2 - 4ac}{2a}}$라는 것은 잘 알려져 있습니다.
하지만 근의 공식은 몇 가지 문제점을 가지고 있는데, 첫째로 5차 이상의 다항식에 대해서는 근의 공식이 존재하지 않는다는 점과, 둘째로 3차 이상의 다항식에 대한 근의 공식으로는 깔끔하게 a + bi의 형태로 해를 타나기 힘들다는 점이 있습니다.
대신 컴퓨터를 이용한 수치적인 방법을 이용하면 주어진 다항식의 정확한 해는 아니더라도, 근사해는 꽤 쉽게 구할 수 있습니다.
Bisection method
가장 쉽게 사용해 볼 수 있는 방법은 Bisection method, 즉 이분법입니다. 다항식은 항상 연속이므로 중간값 정리에 의해, 서로 다른 두 점 a, b에 대해서 $f(a) f(b) < 0$이면 구간 [a, b]에 반드시 $f(x) = 0$이 되는 점 x가 존재하게 됩니다.
이 떄 $c = \frac{a + b}{2}$에 대한 함수값을 구한 뒤, $f(a)f(c) < 0$일 경우 구간을 [a, c]로, $f(b)f(c) < 0$일 경우 구간을 [b, c]로 바꾼 후 반복하면 항상 구간의 길이를 절반으로 줄여나갈 수 있으므로 쉽고 정확하게 해를 찾을 수 있게 됩니다.
하지만 이분법의 최대 단점은 $f(a) f(b) < 0$인 구간 [a, b]를 찾아내야 한다는 점입니다. 여러 x에 대해서 시도해본다면 이러한 구간을 찾아낼 수도 있겠지만, 반드시 찾아낼 수 있다는 보장을 할 수는 없습니다.
Newton’s method
또 다른 잘 알려진 방법으로는 newton’s method가 있습니다. Newton’s method를 간단히 설명하면 특정 점에서 접선을 그어, 해당 접선이 x축과 만난 점을 반복해서 구하다보면 이는 다항식의 해에 수렴하게 되는 방법입니다.
이를 점화식으로 써 보면 다음과 같습니다.
$x_{i+1} = x_i - \frac{f(x_i)}{f’(x_{i})}$
Newton’s method는 Bisection method보다 수렴 속도가 빠른 강력한 방법이지만, $x_0$의 값에 따라 근으로 수렴해가지 않을 수도 있기 때문에 초기값을 잘 골라야 합니다.
Bairstow’s method
이제 bairstow’s method에 대해 알아보겠습니다. Bairstow’s method는 주어진 다항식의 실근 뿐 아니라 허근까지도 알아낼 수 있기 때문에 훨씬 활용도가 높습니다.
Bairstow’s method의 주어진 다항식을 이차식을 인수로 가지도록 하여 계속 쪼개나가는 방법을 사용하였습니다.
주어진 다항식을 $P(x) = a_0 + a_1 x + a_2 x^2 + … + a_n x^n$이라고 합시다.
이 다항식을 $x^2 - rx - s$라는 이차식으로 나누어보면 어떻게 될까요?
$P(x) = Q(x)(x^2 - rx - s) + ax + b$처럼, 일차식이 나머지로 남게 됩니다.
이 때, 만약 $x^2 - rx - s$의 해가 주어진 다항식 $P(x)$의 해와 같다면, $P(x)$는 $x^2 - rx - s$로 정확히 나누어떨어지게 됩니다. 다시 말해, a = 0, b = 0이 됩니다.
따라서, 만약 a, b가 0으로 수렴하도록 r과 s를 잘 조절해나간다면 주어진 다항식 $P(x)$의 근사해를 해로 가지는 이차식을 구할 수 있으며, 몫인 $Q(x)$에 대해서도 이러한 과정을 반복해나가면 다항식의 모든 해를 구할 수 있게 됩니다.
이제 r과 s를 어떻게 조절해나가야 할지를 천천히 살펴보겠습니다. r과 s를 조절하는 과정은 2-dimentional newton method를 사용합니다.
다시 $P(x)$를 $x^2 - rx - s$로 나누는 과정을 한 번 살펴보겠습니다.
이 때 생기는 몫 $Q(x)$와 나머지 $ax + b$의 계수는 어떻게 될까요?
$Q(x) = b_2 + b_3 x + b_4 x^2 + … + b_n x^{n-2}$라고 놓고, 나머지는 $b_0 + b_1(x-r)$이라고 놓을 경우, 계수를 비교해보면 아래와 같은 식을 얻을 수 있습니다.
$b_n = a_n$
$b_{n-1} = a_{n- 1} + r b_n$
$b_{i} = a_{i} + r b_{i+1} + sb_{i+2}$ (i = n - 2 to 0) … 식 1
우리의 목표는 r, s를 변형시켜서 $b_0$와 $b_1$을 0에 가깝게 만드는 것입니다.
즉, $b_0(r, s) \approx 0$, $b_1(r, s) \approx 0$을 만들어야 합니다.
r과 s를 $\Delta r$, $\Delta s$만큼씩 변형시킬 때 $b_0$와 $b_1$은 아래와 같이 변하게 되며, 이 값이 0이 되어야 합니다.
$b_0(r + \Delta r, s + \Delta s) = b_0 + \frac{\partial b_0}{\partial r}\Delta r + \frac{\partial b_0}{\partial s}\Delta s = 0$
$b_1(r + \Delta r, s + \Delta s) = b_1 + \frac{\partial b_1}{\partial r}\Delta r + \frac{\partial b_1}{\partial s}\Delta s = 0$
이제 $\Delta r$, $\Delta s$에 대한 두 개의 방정식을 얻었습니다. 이 연립방정식을 풀기만 하면 $\Delta r$, $\Delta s$를 구할 수 있고, 이에 맞게 r, s를 변형시키면 됩니다.
하지만 아직 한 가지 과정이 더 남아있는데, 바로 연립방정식의 계수를 구하는 과정입니다. 연립방정식의 계수인 $\frac{\partial b_0}{\partial r}$, $\frac{\partial b_0}{\partial s}$, $\frac{\partial b_1}{\partial r}$, $\frac{\partial b_1}{\partial s}$를 구해야 합니다.
이 계수 또한 앞과 유사한 방법을 사용해 구할 수 있는데, 바로 위 식의 $b_i$에 대해 $x^2 - rx - s$로 한 번 더 나누어주는 것입니다.
이 경우 식 1에서 $a_i$를 $b_i$로, $b_i$를 $c_i$로 바꾼 것과 같은 점화식을 가지게 되며 최종적인 식 $c_1 = b_1 + rc_2 + sc_3$, $c_0 = b_0 + rc_1 + sc_2$에서 $\frac{\partial b_0}{\partial r} = c_1$, $\frac{\partial b_0}{\partial s} = \frac{\partial b_1}{\partial r} = c_2$, $\frac{\partial b_1}{\partial s} = c_3$이라는 식을 얻을 수 있습니다.
식이 조금 복잡해졌는데, 최종적으로 정리하면 $c_1 \Delta r + c_2 \Delta s = -b_0$, $c_2 \Delta r + c_3 \Delta s = -b_1$ 두 개의 연립방정식을 풀면 됩니다.
이러한 과정을 반복한다면 최종적으로 알맞은 r과 s를 찾을 수 있게 되며, 주어진 다항식을 이차 이하의 식들의 곱으로 만들 수 있습니다. 이제 각각의 식에 대해 근의 공식을 통해 근을 계산하기만 하면 모든 과정이 완료됩니다.
예제 코드
Bairstow’s method를 C++를 이용하여 구현해 본 코드입니다. 연립방정식을 풀 때 몇 가지 예외처리를 해야 하는데 이는 생략하였습니다.
#include <iostream>
#include <vector>
#include <complex>
#include <cmath>
using namespace std;
vector<complex<double>> bairstow(int n, vector<double> a) {
if (n <= 2) {
if (n == 1) return { {-a[0] / a[1], 0} };
else {
double D = a[1] * a[1] - 4 * a[0] * a[2];
if (D >= 0) {
return { {(-a[1] + sqrt(D)) / (2 * a[2]), 0}, {(-a[1] - sqrt(D)) / (2 * a[2]), 0} };
}
else {
return { {-a[1] / (2 * a[2]), sqrt(-D) / (2 * a[2])}, {-a[1] / (2 * a[2]), -sqrt(-D) / (2 * a[2])} };
}
}
}
double r = -3, s = -3;
vector<double> b(n + 1), c(n + 1);
for (int iter = 0; iter < 1000; ++iter) {
b[n] = a[n];
b[n - 1] = a[n - 1] + r * b[n];
for (int i = n - 2; i >= 0; --i) {
b[i] = a[i] + r * b[i + 1] + s * b[i + 2];
}
c[n] = b[n];
c[n - 1] = b[n - 1] + r * c[n];
for (int i = n - 2; i >= 1; --i) {
c[i] = b[i] + r * c[i + 1] + s * c[i + 2];
}
double dr = (b[0] * c[3] - b[1] * c[2]) / (c[2] * c[2] - c[1] * c[3]);
double ds = (b[1] * c[1] - b[0] * c[2]) / (c[2] * c[2] - c[1] * c[3]);
r += dr;
s += ds;
}
vector<complex<double>> ret(n);
vector<complex<double>> ans1 = bairstow(2, {-s, -r, 1});
vector<complex<double>> ans2 = bairstow(n - 2, vector<double>(b.begin() + 2, b.end()));
for (int i = 0; i < 2; ++i) ret[i] = ans1[i];
for (int i = 2; i < n; ++i) ret[i] = ans2[i - 2];
return ret;
}
int main() {
int n;
cin >> n;
vector<double> a(n + 1);
for (int i = 0; i <= n; ++i) cin >> a[i];
vector<complex<double>> ans = bairstow(n, a);
for (complex<double> u: ans) printf("%lf + %lfi\n", u.real(), u.imag());
}
$f(x) = 6x^5 + 11x^4 - 33x^3 - 33x^2 + 11x + 6$에 대해서, x = -3, -1, $-\frac{1}{3}$, $\frac{1}{2}$, 2라는 5개의 해를 잘 찾아내는 것을 확인할 수 있습니다.
(입력: 5 6 11 -33 -33 11 6)
Reference
다음은 이번 글을 작성하는 데에 참고한 자료들입니다.