一个 O(n^3) 计算矩阵特征多项式的算法

先看一个粗暴的朴素算法

首先是矩阵特征多项式的定义：

$\varphi_A(\lambda) = \det(A - \lambda I)$

最直接的做法就是直接按行展开矩阵行列式并递归计算。这样是一个 $O(n!)$ 的暴力算法，显然是不能接受的。

再看一个多项式复杂度的简单算法

由朴素算法可知，$\varphi_A(\lambda)$ 是一个 $n$ 阶多项式，于是可以带入 $n+1$ 个点进行插值（例如 $\lambda = 0...n$ ）。然后 $O(n^3)$ 地求出一个矩阵的行列式，总体的复杂度就是 $O(n^4)$ 了。

最后是一个优秀的 n^3 算法

先是化简矩阵的部分

由于我本人是一个脏脏的 CS 人，所以我并不打算从 householder-transform 的角度下手 其实是看不懂。

由相似的性质知，若 $B = PAP^{-1}$ 则 $A,B$ 的特征多项式相同。如果能使 $B$ 满足某些性质能更轻松地求出特征多项式，就能借助求解 $B$ 的特征多项式得到 $A$ 的特征多项式。

于是考虑通过初等变换来构造 $P$：

1. 行变化：$E_{i,j} A E_{i,j}^{-1}$ ，互换第 $i,j$ 行、第 $i,j$ 列。
2. 行倍乘：$E_i(k) A E_{i}^{-1}(k)$ ，第 $i$ 行乘以 $k$，第 $i$ 列乘以 $\frac{1}{k}$
3. 行倍加：$E_{i,j}(k) A E_{i,j}^{-1}(k)$，第 $i$ 行加上 $k$ 倍的第 $j$ 行，第 $j$ 列减去 $k$ 倍的第 $i$ 列。

由于每次对第 $i$ 行的操作都会影响到第 $i$ 列，所以考虑将 $A$ 化简为 Upper-Hessenberg 矩阵 上海堡垒矩阵，而非常见的上三角矩阵。

$\begin{matrix} a_{1,1} & a_{1,2} & a_{1,3} & \dots & a_{1,n-1} & a_{1,n} \\ a_{2,1} & a_{2,2} & a_{2,3} & \dots & a_{2,n-1} & a_{2,n} \\ 0 & a_{3,2} & a_{3,3} & \dots & a_{3,n-1} & a_{3,n} \\ 0 & 0 & a_{4,3} & \dots & a_{4,n-1} & a_{4,n} \\ \vdots & \vdots & \vdots & & \vdots & \vdots \\ 0 & 0 & 0 & \dots & a_{n-1,n} & a_{n,n} \end{matrix}$

这样处理第 $i$ 列时，对应的主元在第 $i+1$ 行，就不会因为行倍加操作影响到第 $i$ 列了。

那么如何计算上海森堡矩阵的特征多项式呢

由于 $B$ 和 $A$ 特征多项式一致，所以我们只要算出 $B$ 的特征多项式即可。

令第 $i$ 行到第 $n$ 行，第 $i$ 列到第 $n$ 列相交形成的矩阵的特征多项式为 $D_i(x)$ 。特别的，令 $D_{n+1}(x) = 1$。则按照第一列展开可以得到：

$D_i(x) = (a_{i,i} - x) D_{i+1}(x) + \sum_{j = i + 1}^na_{i, j}(-1)^{j - i}D_{j+1}(x) \prod_{k=i+1}^{j} a_{k,k-1}$

这样直接后往前计算 $D_i(x)$ ，复杂度也是 $O(n^3)$ 的。加上前面化简矩阵，整体的复杂度就是 $O(n^3)$ 了。

一个参考代码：

#include <iostream>
#include <vector>
#include <algorithm>

typedef std::vector<std::vector<double>> matrix;
typedef std::vector<double> polynomial;

void swapLine(matrix &A, int i, int j) {
    std::swap(A[i], A[j]);
}

void swapColumn(matrix &A, int i, int j) {
    for (auto& k : A) {
        std::swap(k[i], k[j]);
    }
}

void addLine(matrix &A, int i, int j, double k) {
    for (int l = 0; l < A[0].size(); l++) {
        A[i][l] += k * A[j][l];
    }
}

void addColumn(matrix &A, int i, int j, double k) {
    for (auto& l : A) {
        l[i] += k * l[j];
    }
}

matrix simplify(matrix A) {
    int n = static_cast<int>(A.size());

    for (int i = 0; i + 1 < n; i++) {
        for (int j = i + 1; j < n; j++) {
            if (A[j][i] != 0) {
                if (j != i + 1) {
                    swapLine(A, i + 1, j);
                    swapColumn(A, i + 1, j);
                }
                break;
            }
        }

        if (A[i + 1][i] == 0) {
            continue;
        }

        for (int k = i + 2; k < n; k++) {
            double t = A[k][i] / A[i + 1][i];
            addLine(A, k, i + 1, -t);
            addColumn(A, i + 1, k, t);
        }
    }

    return A;
}

polynomial calcPolynomial(const matrix& A) {
    int n = static_cast<int>(A.size());
    auto D = std::vector(n + 1, polynomial());

    D[n] = {1};

    for (int i = n - 1; i >= 0; i--) {
        D[i].resize(n - i + 1);

        double k = 1;
        for (int j = i + 1; j < n; j++) {
            k *= A[j][j - 1];

            for (int l = 0; l < D[j + 1].size(); l++) {
                D[i][l] += ((j - i) % 2 == 0 ? 1 : -1) * A[i][j] * k * D[j + 1][l];
            }
        }

        for (int j = 0; j < D[i + 1].size(); j++) {
            D[i][j] += A[i][i] * D[i + 1][j];
            D[i][j + 1] -= D[i + 1][j];
        }
    }

    return D[0];
}

int main() {
    matrix m = {
        {1, 1, 1},
        {1, 2, 3},
        {1, 3, 6}
    };

    m = simplify(m);
    auto p = calcPolynomial(m);

    for (int i = 0; i < p.size(); i++) {
        if (i == 0) {
            std::cout << p[i];
        } else {
            std::cout << " + " << p[i] << "x^" << i;
        }
    }

    std::cout << std::endl;

    return 0;
}