Linear Algebra in Problem Solving (2)

Linear Algebra in Problem Solving (1)

Linear Algebra in Problem Solving (2) (현 포스트)

Linear Algebra in Problem Solving (3)

 

이 포스트에서는 선형대수학 시리즈의 전편에 이어 좀더 고급 알고리즘에 대해 다룬다. 전편에서 만든 코드베이스를 그대로 사용할 예정이므로 전편을 먼저 읽는 것을 권장한다.

특성다항식

\(n\times n\) 행렬 \(A\)의 특성다항식(characteristic polynomial) \(f_A\)는 아래와 같이 정의된다.

\[f_A(x):=\det(xI-A)\]

이 다항식은 최고차항의 계수가 1인 \(n\)차 다항식이며, 아래와 같이 여러 성질을 가지고 있는 중요한 다항식이다.

• \(f_A(x)=0\)의 근은 \(A\)의 고윳값이다.
• \(f_A\)의 \(x^{n-1}\)의 계수는 \(A\)의 trace의 부호를 뒤집은 값이다.
• \(f_A\)의 상수항은 \((-1)^n\det(A)\)이다.
• \(B=P^{-1}AP\)라면, \(f_A=f_B\)이다.
• \(f_A(A)=O\)이다. (Cayley–Hamilton 정리)

마지막의 경우 숫자가 들어가야 할 자리에 행렬이 들어가는 것이 이상하게 생각될 수도 있지만, 행렬의 거듭제곱과 상수 곱셈도 잘 정의되므로 그런 맥락으로 이해하면 된다.

 

이 특성다항식은 어떻게 계산할 수 있을까? 전편에서 알아본 행렬식 구하는 알고리즘으로 구하게 되면, 행렬의 원소가 함수이기 때문에 더이상 덧셈과 곱셈이 상수 시간에 동작하지 않으며, 설상가상으로 나눗셈도 까다로워진다. 따라서 다른 방법을 생각해야 한다.

 

핵심은 위 특성다항식의 성질 중 네 번째, 즉 \(B=P^{-1}AP\)이면 \(f_A=f_B\)라는 사실을 이용하는 것이다. 저런 \(P\)가 존재하면 \(A\)와 \(B\)를 닮은(similar) 행렬이라고 하는데, \(A\)와 닮음이면서 형태가 간단해서 특성다항식을 계산하기 쉬운 \(B\)를 찾을 수 있다면, 우리는 \(A\) 대신 \(B\)의 특성다항식을 계산하기만 하면 된다.

 

그럼 특성다항식을 계산하기 쉬운 행렬은 어떤게 있을까? 대각행렬과 삼각행렬은 특성다항식을 계산하기 매우 쉽다. 하지만 모든 행렬이 대각화나 삼각화가 가능하지는 않고, 설사 가능하다고 해도 계산하는 것이 쉽지 않다. 따라서 대안으로 사용하는 것이 Hessenberg matrix이다. 정확히는 우리는 upper Hessenberg matrix를 사용할 것인데, upper triangular matrix가 주대각선 아래로 모두 0인 행렬이라면, upper Hessenberg matrix는 마찬가지로 주대각선 아래로 대부분 0인 행렬이지만, 예외적으로 주대각선 바로 아래까지는 0이 아니어도 되는 행렬이다. 예를 들어 다음과 같은 행렬이 upper Hessenberg matrix이다.

\[\begin{bmatrix}1&5&5&0&-2\\ -1&2&3&7&1\\ 0&4&4&-2&0\\ 0&0&3&2&2\\ 0&0&0&7&-1\end{bmatrix}\]

수학적으로 표현하면, \(B\)의 \(i\)번째 행 \(j\)번째 열에 오는 원소를 \(b_{ij}\)라고 할 때, \(i\geq j+2\)일 때 \(b_{ij}=0\)이면 \(B\)를 upper Hessenberg matrix라고 한다.

 

주어진 행렬과 닮음인 upper Hessenberg matrix를 찾는 방법을 논하기 전에, upper Hessenberg matrix의 특성다항식을 구하는 방법을 먼저 알아보자. 어떤 행렬 \(M\)의 첫 \(k\)개의 행과 열로 이루어진 부분행렬을 leading principal minor라고 하는데, 이 행렬을 \([M]_k\)라고 쓰자. \(M\)이 upper Hessenberg matrix라면 모든 leading principal minor도 upper Hessenberg이기 때문에, \(\det[M]_k\)를 \(k\)번째 열에 cofactor expansion을 적용하여 재귀적으로 구할 수 있다. 식으로 적으면 아래와 같다.

\[\det[M]_k=\sum_{i=0}^{k-1}(-1)^{k-1-i}\det[M]_i\cdot m_{i+1,k}\prod_{j=i+2}^{k}m_{j,j-1}\]

여기서 \(m_{ij}\)는 \(M\)의 \(i\)번째 행 \(j\)번째 열에 있는 원소이고, \(\det[M]_0\)은 1로 정의한다.

 

이렇게 계산하면 동적 계획법(dynamic programming)으로 \(\det(M)=\det[M]_n\)을 계산하는 데 \(O(n^2)\)의 시간이면 충분하기 때문에, 일반적인 행렬의 행렬식을 계산하는 것보다 훨씬 빠르다. 이제 \(B\)가 upper Hessenberg matrix이면 \(xI-B\)도 upper Hessenberg이기 때문에, 위 점화식에서 \(M=xI-B\)를 대입해 보자. 이 행렬은 원소 중에 일차식이 있기는 하지만 주대각선에만 일차식이 있고 나머지는 모두 상수항이다. 따라서 점화식에서 \(\det[M]_i\) 뒤에 오는 부분은, \(i=k-1\)일 때 \(m_{i+1,k}\)만 일차식이고 나머지는 모조리 상수이다! 따라서 \(\det[M]_k\)을 구할 때 \(k\)차 이하의 다항식 \(k\)개를 더하는 것이 시간을 가장 오래 차지하는 부분이고, 총 \(O(n^3)\)의 시간이면 충분하다는 것을 알 수 있다.

 

정리하자면, (upper) Hessenberg matrix는 \(O(n^3)\) 시간에 특성다항식을 계산 가능하다.

 

이제 남은 것은 임의의 행렬 \(A\)와 닮음인 upper Hessenberg matrix \(B\)를 찾는 것이다. 이를 위해서는 elementary row operation을 응용하면 되는데, 아쉽게도 elementary row operation이 적용된 행렬은 원래 행렬과 닮음이 아닐 수 있다. 하지만 elementary row operation을 적용한 후에 비슷한 elementary column operation을 적용하여 원래 행렬과 닮음이 되게 할 수 있다! 각 연산별로 어떻게 하는지 알아보자.

 

먼저 가장 간단한 row switching 연산의 경우, \(i\)번째 행과 \(j\)번째 행을 바꾸는 것은 아래와 같이 permutation matrix를 왼쪽에 곱한 것과 같다.

\[\begin{bmatrix}1&0&0&0&0\\0&0&0&1&0\\0&0&1&0&0\\0&1&0&0&0\\0&0&0&0&1\end{bmatrix}A\]

예시는 2번째 행과 4번째 행을 바꾸는 경우이다. 이후 \(A\)와 닮음으로 만들어주기 위해서는 해당 permutation matrix의 역행렬을 오른쪽에 곱해줘야 한다. 두 위치만 바꾸는 연산은 역연산이 자기 자신이므로, 같은 행렬을 오른쪽에 곱해주면 된다.

\[\begin{bmatrix}1&0&0&0&0\\0&0&0&1&0\\0&0&1&0&0\\0&1&0&0&0\\0&0&0&0&1\end{bmatrix}A\begin{bmatrix}1&0&0&0&0\\0&0&0&1&0\\0&0&1&0&0\\0&1&0&0&0\\0&0&0&0&1\end{bmatrix}\]

오른쪽에 이 행렬을 곱해주는 것은, \(i\)번째 열과 \(j\)번째 열을 바꾸는 것과 효과가 같다. 따라서 row switching 연산은, 같은 인덱스를 바꾸는 column switching 연산이 뒤따라 오면 원래 행렬과 닮음이 된다.

 

다음으로 row addition 연산의 경우, \(i\)번째 행에 \(j\)번째 행의 \(c\)배를 더하는 것은 아래와 같이 곱하는 것과 같다.

\[\begin{bmatrix}1&0&0&0&0\\0&1&0&0&0\\0&0&1&0&0\\0&0&0&1&0\\0&-2&0&0&1\end{bmatrix}A\]

예시는 5번째 행에 2번째 행의 -2배를 더하는 연산이다. 이 행렬의 역행렬은 \(c\)의 부호를 뒤집는 것으로, 닮음으로 만들기 위해 아래와 같이 오른쪽에 곱해주면 된다.

\[\begin{bmatrix}1&0&0&0&0\\0&1&0&0&0\\0&0&1&0&0\\0&0&0&1&0\\0&-2&0&0&1\end{bmatrix}A\begin{bmatrix}1&0&0&0&0\\0&1&0&0&0\\0&0&1&0&0\\0&0&0&1&0\\0&2&0&0&1\end{bmatrix}\]

이것은 \(j\)번째 열에 \(i\)번째 열의 \(-c\)배를 더하는 것과 같은 효과이다. 따라서 row addition 연산은, 부호를 바꾸고 두 인덱스의 역할을 바꾼 column addition 연산이 뒤따라 오면 원래 행렬과 닮음이 된다.

 

row multiplication이 빠졌는데, 쉬울 뿐더러 실제로 활용하지 않을 예정이므로 생략한다.

 

위 사실을 이용하면, 가우스 소거법과 흡사한 방식으로 \(A\)와 닮음인 upper Hessenberg matrix를 만들 수 있다.

1. 1열부터 \(n-2\)번째 열까지 순서대로 처리한다.
2. 현재 \(j\)번째 열에 있다면, \((j+1)\)번째 행부터 내려가면서 0이 아닌 원소의 위치를 찾는다.
3. 그러한 원소가 없다면 다음 열로 넘어가서 2번 과정으로 돌아간다.
4. 그러한 원소가 발견되면 row switching으로 \((j+1)\)번째 행으로 가져오고, 대응되는 column switching을 적용한다.
5. \((j+1)\)번째 행에 적절한 상수를 곱한 후 더해서 \((j+2)\)번째 행부터 모두 0이 되도록 row addition 연산을 적용하고, 각 연산 이후에 대응되는 column addition 연산도 수행한다.

이 과정을 잘 살펴보면, 닮음을 유지하기 위한 elementary column operation이 현재 작업중인 열에 영향을 주지 않는다는 것을 알 수 있다. 따라서 모든 과정이 끝나면 원래 행렬과 닮음인 upper Hessenberg matrix가 완성된다.

 

이 과정의 시간복잡도는 가우스 소거법과 동일한 \(O(n^3)\)이고, Hessenberg matrix의 특성다항식 알고리즘을 적용하면 전체 \(O(n^3)\) 시간에 임의의 행렬의 특성다항식을 구할 수 있게 된다.

 

이제 이 과정을 구현해보자. 먼저 아래와 같이 elementary column operation을 추가한다. 이 코드에서는 column switching 연산을 column index 배열을 이용하는 트릭을 사용하여 상수 시간에 구현하였다. Elementary row operation 부분은 동일하므로 생략하였다.

template<class T>
class matrix {
public:
    int n, m;
    vector<int> cidx;
    vector<vector<T>> entry;

    matrix(int n_, int m_) : n(n_), m(m_), cidx(m_), entry(n_, vector<T>(m_)) {
        iota(cidx.begin(), cidx.end(), 0);
    }
    matrix(int n_) : matrix(n_, n_) {}

    T& operator[](const pair<int, int>& p) {
        return entry[p.first][cidx[p.second]];
    }

    const T& operator[](const pair<int, int>& p) const {
        return entry[p.first][cidx[p.second]];
    }

    /* ... */

    void cswap(int i, int j) {
        swap(cidx[i], cidx[j]);
    }

    void cmul(int j, T x) {
        for (int i = 0; i < n; i++) {
            entry[i][cidx[j]] *= x;
        }
    }

    void cadd(int j1, int j2, T x) {
        for (int i = 0; i < n; i++) {
            entry[i][cidx[j1]] += entry[i][cidx[j2]] * x;
        }
    }
};

이후 이 연산을 이용해서 특성다항식을 구하면 된다. Column multiplication 연산은 사용하지 않는다.

template<class T>
vector<T> char_poly(matrix<T> &A) {
    int n = A.n;
    assert(n == A.m);

    /* Transform to upper Hessenberg matrix */
    for (int j = 0; j < n - 2; j++) {
        for (int i = j + 1; i < n; i++) {
            if (A[{i, j}] != 0) {
                A.rswap(i, j + 1);
                A.cswap(i, j + 1);
                break;
            }
        }
        if (A[{j + 1, j}] != 0) {
            for (int i = j + 2; i < n; i++) {
                T x = -A[{i, j}] / A[{j + 1, j}];
                A.radd(i, j + 1, x);
                A.cadd(j + 1, i, -x);
            }
        }
    }

    /* Compute the char poly */
    vector<vector<T>> dp(n + 1);
    dp[0] = {1};
    for (int i = 1; i <= n; i++) {
        dp[i].resize(i + 1);
        for (int j = 1; j <= i; j++) {
            dp[i][j] = dp[i - 1][j - 1];
        }
        for (int j = 0; j < i; j++) {
            dp[i][j] -= A[{i - 1, i - 1}] * dp[i - 1][j];
        }
        T prod = 1;
        for (int j = i - 2; j >= 0; j--) {
            prod *= A[{j + 1, j}];
            for (int jj = 0; jj <= j; jj++) {
                dp[i][jj] -= A[{j, i - 1}] * prod * dp[j][jj];
            }
        }
    }
    return dp[n];
}

가우스 소거법 구현과 마찬가지로 \(\mathbb{F}_2\)에서 bitset을 사용할 수 있지만, 열 연산의 존재로 인해 행렬을 \(n^2\) 크기의 bitset으로 구현해야 해서 구현은 좀더 까다롭다. 이 경우 시간복잡도는 \(O(n^3/w)\).

최소다항식

앞서 짧게 언급한 Cayley–Hamilton 정리에 의하면 \(f_A(A)=0\)이 된다. 최소다항식 \(p_A(x)\)는 \(p_A(A)=O\)을 만족하는 다항식 중 가장 차수가 낮으면서 최고차항의 계수가 1인 다항식이다. 수학적으로 이러한 다항식이 유일함을 증명할 수 있으며, \(g(A)=O\)을 만족하는 모든 다항식 \(g(x)\)가 \(p_A(x)\)로 나누어 떨어진다는 사실도 증명할 수 있다. 따름정리로 \(f_A(x)\) 역시 \(p_A(x)\)로 나누어떨어진다.

 

PS에서 최소다항식을 사용해야하는 경우는 매우 드물지만, 최소다항식을 구하는 알고리즘 역시 가우스 소거법의 좋은 적용 예이므로 소개하고자 한다. 핵심 아이디어는 행렬을 \(n^2\)차원 벡터로 취급하는 것이다. 그럼 아래 나열된 행렬의 선형 결합이 0이 되는 계수 조합 중 0이 아닌 계수가 가능한 한 앞쪽에 있게 하고 싶은 것이다.

\[I, A, A^2, \cdots, A^n\]

이제 아래와 같이 열벡터가 위 행렬인 \(n^2\times (n+1)\) 행렬을 만들고 RREF 변환 과정(단순 가우스 소거법이 아니다!)을 진행하다가, 처음으로 leading one을 찾을 수 없는 열이 나오자마자 멈추면 된다. 행렬과 구분하기 위해 행렬을 \(n^2\)차원 열벡터로 편 것을 굵게 표시하였다.

\[\left[\begin{array}{c|c|c|c|c}\mathbf{I}&\mathbf{A}&\mathbf{A^2}&\cdots&\mathbf{A^n}\end{array}\right]\]

그러면 그러한 열이 나왔을 때 최소다항식의 계수는 어떻게 될까? 그 열벡터에 적혀있다! 정확히는 \(k\)번째 열벡터에서 멈췄다면 해당 열벡터의 첫번째 원소부터 \((k-1)\)번째 원소까지에 최소다항식의 상수항부터 \(x^{k-2}\)의 계수까지가 부호가 뒤집힌 채로 적혀있다. 최고차항의 계수는 당연히 1이다.

template<class T>
vector<T> min_poly(const matrix<T>& A) {
    int n = A.n;
    assert(n == A.m);
    vector<matrix<T>> mv(n + 1, matrix<T>(n));
    for (int i = 0; i < n; i++) {
        mv[0][{i, i}] = 1;
    }
    for (int i = 1; i <= n; i++) {
        mv[i] = matmul(mv[i - 1], A);
    }
    matrix<T> M(n * n, n + 1);
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            for (int k = 0; k <= n; k++) {
                M[{i * n + j, k}] = mv[k][{i, j}];
            }
        }
    }
    rref(M);
    for (int d = 1; d <= n; d++) {
        if (M[{d, d}] == 0) {
            vector<T> res(d + 1);
            for (int i = 0; i < d; i++) {
                res[i] = -M[{i, d}];
            }
            res[d] = 1;
            return res;
        }
    }
    assert(0);
    return vector<T>();
}

시간복잡도는 \(O(n^4)\)이며, \(\mathbb{F}_2\)에서 bitset을 활용할 경우 \(O(n^4/w)\)로 줄일 수 있다.

\(\det(xA+B)\)

\(n\times n\) 행렬 \(A\)와 \(B\)가 있을 때, 다항식 \(\det(xA+B)\)를 구하는 문제를 생각해 보자. 이 형태는 특히 조합론 문제를 선형대수학을 이용해 푸는 예시에서 종종 볼 수 있다. (이러한 예시는 다음 포스트에서 좀더 자세히 다룰 예정이다.) 만약 \(A\)가 역행렬이 존재한다면 \(-A^{-1}B\)의 특성다항식을 구하는 문제와 동일하므로 특성다항식의 확장이라고 볼 수도 있다. 하지만 \(A\)가 역행렬이 존재하지 않을 수도 있기 때문에 특성다항식을 구할 때 사용한 방식을 적용하려면 추가적인 처리가 필요하다. 이러한 처리를 하는 방법에는 여러 가지가 있지만, 그중 한 가지를 소개한다.

 

목표는 \(xA+B\)에 elementary row operation과 elementary column operation을 적용하여 아래와 같은 형태로 만드는 것이다.

\[\left[\begin{array}{@{}c|c@{}}\begin{matrix}\\xI_k+C\\ \quad \end{matrix} & O \\ \hline  D &  L\end{array}\right]\]

여기서 \(I_k\)는 \(k\times k\) 단위행렬, \(C\)는 상수로만 이루어진 행렬, \(O\)는 영행렬, 그리고 \(L\)은 주대각선이 모두 1인 \((n-k)\times (n-k)\) lower triangular matrix이다. 이런 행렬의 행렬식은

\[\det(xI+C)\cdot\det(L)=\det(xI+C)\]

이므로 \(-C\)의 특성다항식을 구하면 된다. 즉 우리는 기본 연산만으로 \(xA+B\)를 이런 형태로 바꾸는 방법을 알면 된다.

 

이 변환 과정을 재귀적으로 생각해보면 어렵지 않게 방법을 찾을 수 있다. 첫 \(\alpha\)개의 열과 마지막 \(\beta\)개의 열이 우리의 의도대로 이미 변형되었다고 가정하고, 그 사이의 열에서 작업을 하는 것이 핵심이다.

1. 처음에 \(\alpha=\beta=0\)이다.
2. \(\alpha+\beta=n\)이 되면 처리가 완료되었으므로 종료한다.
3. \((\alpha+1)\)번째 열을 본다. 첫번 째 행부터 \(\alpha\)번째 행까지에 일차항이 있으면 첫번째 열부터 \(\alpha\)번째 열까지를 이용해 column addition 연산으로 일차항을 모두 없애준다.
4. \((\alpha+1)\)번째 열에서 \((\alpha+1)\)번째 행부터 \((n-\beta)\)번째 행까지 순서대로 보면서 일차항이 있는지 찾는다.
5. 일차항을 발견하면 elementary row operation을 적절히 이용해 \((\alpha+1\)번째 행에만 일차항 \(x\)가 남아있게 한다. \(\alpha\)에 1을 더하고 2번으로 돌아간다.
6. 일차항을 발견하지 못하면 \((n-\beta)\)번째 행까지는 상수만 있다는 뜻이다. 해당 열을 column switching으로 \((n-\beta)\)번째 열로 옮겨준다.
7. 이제 \((n-\beta)\)번째 열을 본다. \((n-\beta)\)번째 행부터 역순으로 보면서 0이 아닌 원소가 있는지 찾는다.
8. 없으면 행렬식이 0이라는 뜻이므로 바로 종료해도 된다.
9. 0이 아닌 원소를 찾았고 행의 번호가 \((\alpha+1)\) 이상일 경우 row switching으로 \((n-\beta)\)번째 행으로 옮긴 후 row multplication으로 해당 원소를 1로 만들고 row addition으로 그 이전 행까지를 모두 0으로 만든다. 이후 \(\beta\)에 1을 더하고 2번으로 돌아간다.
10. 0이 아닌 원소를 찾았고 행의 번호가 \(\alpha\) 이하인 경우 row switching으로 \(\alpha\)번째 행으로 옮겨오고 같은 인덱스로 column switching까지 진행한다.
11. 이제 row switching으로 \((n-\beta)\)번째 행까지 옮기고 9번과 유사하게 처리한다. 이후 \(\beta\)에 1을 더하지만, \(\alpha\)번째 열의 정합성이 깨졌을 수 있으므로 \(\alpha\)에서 1을 뺀 후 2번으로 돌아간다.

과정이 길어 복잡하게 느껴지지만, 직접 행렬을 그려가면서 따라가면 직관적으로 이해할 수 있다. 10번 과정이 주목할만한데, column switching을 수행하는 이유는 첫 열부터 \(\alpha\)번째 열까지의 정합성을 깨뜨리는 것이 불가피하기 때문에 그것을 마지막인 \(\alpha\)번째 열로 설정하기 위해서이다.

template<class T>
vector<T> linear_det(matrix<T>& A, matrix<T>& B) {
    int n = A.n;
    assert(n == A.m && n == B.n && n == B.m);
    T det = 1;

    /* ... */

    int a = 0, b = 0;
    while (a + b < n) {
        for (int i = 0; i < a; i++) {
            cadd(a, i, -A[{i, a}]);
        }
        for (int i = a; i < n - b; i++) {
            if (A[{i, a}] != 0) {
                rswap(a, i);
                break;
            }
        }
        if (A[{a, a}] != 0) {
            rmul(a, 1 / A[{a, a}]);
            for (int i = a + 1; i < n; i++) {
                radd(i, a, -A[{i, a}]);
            }
            a++;
        } else {
            int r = n - 1 - b;
            cswap(a, r);
            int pos = -1;
            for (int i = r; i >= 0; i--) {
                if (B[{i, r}] != 0) {
                    pos = i;
                    break;
                }
            }
            if (pos == -1) {
                return vector<T>();
            } else {
                if (pos < a) {
                    rswap(pos, a - 1);
                    cswap(pos, a - 1);
                    rswap(a - 1, r);
                    a--;
                } else {
                    rswap(pos, r);
                }
                rmul(r, 1 / B[{r, r}]);
                for (int i = 0; i < r; i++) {
                    radd(i, r, -B[{i, r}]);
                }
                b++;
            }
        }
    }
    matrix<T> C(a);
    for (int i = 0; i < a; i++) {
        for (int j = 0; j < a; j++) {
            C[{i, j}] = -B[{i, j}];
        }
    }
    auto res = char_poly(C);
    for (int i = 0; i <= a; i++) {
        res[i] *= det;
    }
    return res;
}

주석 처리된 부분에는 원래 코드에 쓰이는 행 연산과 열 연산 함수를 \(A\)와 \(B\)에 모두 적용하고 그에 맞게 행렬식 변수를 조정해주는 보조 함수가 있어야 하나, 가독성을 위해 생략하였다.

 

이제 시간복잡도를 확인해보자. 메인 루프 내의 모든 작업은 \(O(n^2)\)에 이루어지고, 루프를 돌 때마다 \(\alpha+2\beta\)의 값이 항상 증가하므로 루프를 도는 횟수는 \(O(n)\)이다. 따라서 전체 시간복잡도는 \(O(n^3)\)으로 특성다항식을 구하는 시간복잡도와 같다. \(\mathbb{F}_2\)에서 bitset을 이용할 경우 마찬가지로 \(O(n^3/w)\).

종결식

종결식(resultant)은 선형대수학의 범주에 속하는 것은 아니라고 볼 수도 있으나, 고난이도 선형대수학 문제에서 종종 등장하므로 간략하게 짚고 넘어가려고 한다.

 

종결식은 두 다항식 \(A(x)=a_0+a_1x+\cdots+a_dx^d\)와 \(B(x)=b_0+b_1x+\cdots+b_ex^e\)에 대해서 정의되며, 아래와 같다(위키피디아와 계수 인덱스가 다르니 주의).

\[\operatorname{res}(A, B):=\det\begin{bmatrix} a_d      & 0           & \cdots & 0          & b_e        & 0              & \cdots & 0       \\ a_{d-1}    & a_d       & \cdots & 0           & b_{e-1}     & b_e           & \cdots & 0  \\ a_{d-2}    & a_{d-1}     & \ddots & 0           & b_{e-2}     & b_{e-1}         & \ddots & 0 \\ \vdots  &\vdots   & \ddots & a_d        & \vdots   &\vdots       & \ddots & b_e  \\ a_0       & a_1 & \cdots & \vdots   & b_0       & b_1     & \cdots & \vdots\\ 0          & a_0       & \ddots &  \vdots  & 0          & b_0          & \ddots &  \vdots  \\ \vdots  & \vdots   & \ddots & a_1  & \vdots  & \vdots      & \ddots & b_1   \\ 0          & 0          & \cdots  & a_0       & 0           & 0              & \cdots & b_0   \end{bmatrix}\]

만일 두 다항식이 일차식의 곱으로 완전히 인수분해되어 \(A(x)=a_d(x-\alpha_1)\cdots(x-\alpha_d)\)와 \(B(x)=b_e(x-\beta_1)\cdots(x-\beta_e)\)와 같이 표현된다면, 종결식은 아래와 같은 식으로도 구할 수 있다.

\[\operatorname{res}(A,B)=a_d^eb_e^d\prod_{i,j}(\alpha_i-\beta_j)\]

이는 근과 계수와의 관계에 의해 아래와 같이 쓸 수도 있다.

\[\operatorname{res}(A,B)=a_d^e\prod_{i=1}^dB(\alpha_i)=(-1)^{de}b_e^d\prod_{j=1}^eA(\beta_j)\]

현대대수학을 조금 알고 있는 독자라면 행렬식을 이용한 정의는 임의의 가환환(commutative ring)에서 정의되지만, 근을 이용한 정의는 그 환을 포함하는 대수적으로 닫힌 체(algebraically closed field)가 있어야 하므로 integral domain에서만 정의된다는 사실을 알고 있으면 좋다.

 

위 식을 이용해 아래와 같은 여러 성질을 쉽게 얻어낼 수 있다.

• \(\operatorname{res}(A,B)=(-1)^{de}\operatorname{res}(B,A)\)
• \(\operatorname{res}(AB,C)=\operatorname{res}(A,C)\operatorname{res}(B,C)\)
• \(A(x)=0\)와 \(B(x)=0\)가 공통근을 가질 필요충분조건은 \(\operatorname{res}(A,B)=0\)이다.
• \(\operatorname{res}(aA, bB)=a^eb^d\operatorname{res}(A,B)\)
• 차수가 \(r<e\)인 다항식 \(R(x)\)에 대해 \(B=AQ+R\)이라면, \(\operatorname{res}(A,B)=a_d^{e-r}\operatorname{res}(A,R)\)

종결식을 계산할 때 가장 중요한 성질은 마지막 성질로, 유클리드 알고리즘으로 두 다항식의 최대공약수를 구하는 과정에서 두 다항식의 종결식도 쉽게 구할 수 있다는 것이다. 특히 최고차항의 계수가 1이면 \(a_d^{e-r}=1\)이므로 계산하기 더욱 쉽다.

 

두 다항식이 어떤 형태로 주어지느냐에 따라 코드와 시간복잡도가 다양하게 나올 수 있으므로, 종결식을 구하는 코드나 시간복잡도 분석은 제공하지 않는다. 하지만 이 글을 읽는 독자라면 쉽게 구현할 수 있을 것이다.

마치며

이 포스트에서는 조금 더 고급 알고리즘인 특성다항식, 최소다항식 등을 구하는 알고리즘에 대해 알아보았다. 이 정도 알고리즘이 필요한 문제는 BOJ와 solved.ac에서 상당히 고난이도 문제로 취급된다. 좀더 다양한 문제를 풀기 위해서는 조합론이나 그래프 이론 등에서 선형대수학이 응용되는 예시를 함께 알면 좋은데, 이는 다음 포스트에서 다룰 예정이다.

연습문제

난이도 순이 아니다.

BOJ 18163: Binary Matrix

BOJ 18542: Permutant

BOJ 19562: Matrix and Queries

BOJ 25509: Lord of the Characteristic Polynomials (1) (주어진 대수 구조가 integral domain이 아닐 수 있으므로 주의가 필요하다.)