Wiedemann Algorithm

Purpose: Solve a sparse linear system, compute a determinant, or find the minimal polynomial of a sparse matrix over a finite field in $O(n\cdot \mathrm{nnz})$ operations — where $\mathrm{nnz}$ is the number of nonzeros — using only matrix-vector products. No fill-in, no $O(n^3)$.

The Idea

Gaussian elimination on a sparse matrix destroys sparsity: zeros fill in, and you end up doing $O(n^3)$ work on a dense matrix that no longer fits in memory. Wiedemann never modifies the matrix at all. It only ever computes $Av$.

1. Pick random vectors $u,v$.
2. Compute the Krylov sequence of scalars
   $s_k=u^T{A}^{k}v,k=0,1,\dots ,2n-1.$
   Each step is one matrix-vector product: $O(\mathrm{nnz})$ each, $O(n\cdot \mathrm{nnz})$ total.
3. Run Berlekamp-Massey on $(s_k)$ to find the shortest linear recurrence it satisfies. With high probability this recurrence is the minimal polynomial $m(\lambda )$ of $A$.
4. From $m(\lambda )={\lambda }^d+\cdots +m_1\lambda +m_0$:
   • Solve $Ax=b$: if $m_0\ne 0$, then $A^{-1}=-\frac{1}{m_0}(A^{d-1}+\cdots +m_{1}I)$, so $x$ is a short linear combination of Krylov vectors.
   • Determinant: if the minimal polynomial equals the characteristic polynomial (generic case), $\det A=(-1{)}^n{m}_0$; otherwise randomise with a diagonal preconditioner and retry.

Complexity

• Time: $O(n\cdot \mathrm{nnz}+n^2)$ field operations
• Space: $O(n+\mathrm{nnz})$ — the matrix is never modified
• Randomized (Las Vegas): the answer is verifiable, so failure only costs a retry

Compare Gaussian elimination: $O(n^3)$ time and $O(n^2)$ space after fill-in. For $n=10^6$ with $\mathrm{nnz}=10^7$, Wiedemann is feasible and elimination is not.

Why the random vectors work

$u^T{A}^{k}v$ satisfies the linear recurrence given by the minimal polynomial of $A$, for any $u,v$. The risk is that the sequence satisfies a shorter recurrence by accident — which happens only if $u$ or $v$ is unluckily aligned with an invariant subspace. Over a field of size $q$ that has probability $O(n/q)$, so for a large prime modulus one try almost always suffices.

Where it is used

• Integer factorisation. The quadratic sieve and the number field sieve both end with a huge sparse ${\mathbb{F}}_2$ linear algebra problem — millions of rows, a few dozen nonzeros per row. Block Wiedemann and block Lanczos are the algorithms that solve it; this step is a substantial fraction of every factorisation record.
• Discrete logarithms by index calculus — the same sparse linear algebra.
• Sparse determinants and ranks in computer algebra systems.

Block Wiedemann

Coppersmith’s block variant uses $u,v$ as blocks of $b$ vectors rather than single vectors, so the matrix-vector products become matrix-matrix products with better cache behaviour, and — crucially — the work parallelises across machines. This is the standard version in factorisation software.

Variants / Use Cases

• Berlekamp-Massey — the subroutine that finds the recurrence
• Lanczos — the alternative sparse method, better suited to symmetric systems
• Kitamasa / Bostan-Mori — evaluate a linear recurrence quickly once you have it
• Rank and Linear Systems — the topic page