Faddeev-LeVerrier Algorithm

Purpose: Compute the characteristic polynomial $\det (\lambda I-A)$ of an $n\times n$ matrix in $O(n^4)$ using only matrix multiplications and traces — no symbolic determinant expansion, no divisions except by small integers.

The Recurrence

Let $p(\lambda )={\lambda }^n+c_{n-1}{\lambda }^{n-1}+\cdots +c_0$. Set $M_0=I$ and iterate for $k=\mathrm{1..}n$:

$$c_{n-k}=-\frac{1}{k}\mathrm{tr}(A{M}_{k-1}),M_k=A{M}_{k-1}+c_{n-k}I.$$

At the end $M_n=0$, which is a free correctness check (and is exactly the Cayley-Hamilton theorem).

Code

// characteristic polynomial coefficients, highest degree first: [1, c_{n-1}, ..., c_0]
vector<double> faddeevLeVerrier(vector<vector<double>> A) {
    int n = A.size();
    vector<double> c(n + 1, 0);
    c[0] = 1;
    vector<vector<double>> M(n, vector<double>(n, 0));
    for (int i = 0; i < n; i++) M[i][i] = 1;               // M = I
 
    for (int k = 1; k <= n; k++) {
        auto AM = matmul(A, M);
        double tr = 0;
        for (int i = 0; i < n; i++) tr += AM[i][i];
        c[k] = -tr / k;
        M = AM;
        for (int i = 0; i < n; i++) M[i][i] += c[k];        // M = A*M + c_k I
    }
    return c;
}

Complexity

• Time: $O(n^4)$ — $n$ iterations, each one $O(n^3)$ matrix multiplication
• With fast matrix multiplication: $O(n^{\omega +1})$
• Space: $O(n^2)$

Faster alternatives exist: reduction to Hessenberg form followed by a recurrence gives $O(n^3)$, and is what numerical libraries use.

Why it works

The identity behind it is Newton’s relation between power sums and elementary symmetric polynomials. The traces $\mathrm{tr}(A^k)$ are the power sums of the eigenvalues, and the coefficients $c_i$ are (up to sign) the elementary symmetric polynomials of the eigenvalues. Newton’s identities convert between them, and the $M_k$ recurrence is a compact way to accumulate exactly those quantities.

Division by $k$

Over a field of characteristic $p\le n$, dividing by $k$ can fail. In particular this algorithm breaks modulo a small prime. Use the Hessenberg method or interpolation instead when working mod $p$ with $p\le n$.

Free byproducts

• Matrix inverse. From Cayley-Hamilton, $A^{-1}=-\frac{1}{c_0}{M}_{n-1}$ (when $c_0\ne 0$).
• Cayley-Hamilton check. $M_n=0$ verifies the computation.
• Determinant. $\det A=(-1{)}^n{c}_0$.
• Trace of every power $\mathrm{tr}(A^k)$ appears along the way — useful for counting closed walks of each length in a graph.

Contest use: linear recurrences

The characteristic polynomial is the bridge between matrix exponentiation and Kitamasa/Bostan-Mori. Given a transition matrix $A$, its characteristic polynomial gives a length-$n$ linear recurrence for the sequence, after which $A^{k}v$ can be computed in $O(n^2\log k)$ or $O(n\log n\log k)$ instead of $O(n^3\log k)$. For a $100\times 100$ transition matrix and $k=10^{18}$, that is the difference between too slow and comfortable.

Variants / Use Cases

• Berlekamp-Massey — recover the recurrence from the sequence directly, usually easier than computing the characteristic polynomial
• Wiedemann — the sparse-matrix analogue, $O(n\cdot \mathrm{nnz})$
• Determinant and Matrices — the topic pages
• Bareiss — for exact integer determinants