v. Sparse and Structured Matrices

When the matrix is sparse or has structure, never run $O(n^3)$ elimination: black-box methods (Wiedemann, Lanczos) only use products $Av$, and structured matrices (circulant, Vandermonde, tridiagonal, rank-one updates) have closed forms or linear-time solvers. mint from Modular Arithmetic.

Matrix	Task	Method	Cost
sparse, finite field	minimal polynomial, solve, determinant, rank	Wiedemann (Krylov sequence + Berlekamp–Massey)	$O(n\cdot \mathrm{nnz}+n^2)$, Las Vegas
sparse symmetric, reals	extreme eigenvalues, systems	Lanczos (three-term recurrence), conjugate gradient	$O(k\cdot \mathrm{nnz})$ for $k$ steps
circulant $M_{ij}=c_{(j-i)\mathrm{mod}n}$	determinant, eigenvalues $C({\omega }^k)$	resultant $\mathrm{Res}(x^n-1,C)$; multiply via cyclic convolution	$O(n^2)$ (or $O(n\log n)$ with FFT over suitable fields)
“permutant” (row $i+1$ = row $i$ permuted by $p$)	determinant	0 unless $p$ is one $n$-cycle, else a permuted circulant	$O(n^2)$
$xI+a{b}^T$	determinant	$x^n+x^{n-1}{b}^{T}a$	$O(n)$
Vandermonde $x_i^j$	determinant, inverse	${\prod }_{i<j}(x_j-x_i)$; columns of the inverse = Lagrange basis polynomials	$O(n^2)$
tridiagonal	solve	Thomas algorithm (elimination without pivoting)	$O(n)$
$A+u{v}^T$ given $A^{-1}$	inverse update	Sherman–Morrison	$O(n^2)$

Wiedemann

Elimination on a sparse matrix fills it in ($O(n^3)$ time, $O(n^2)$ memory); Wiedemann never modifies the matrix. Pick random $u,v$, compute the scalars $s_k=u^T{A}^{k}v$ for $k<2n$ (one sparse product each), and run Berlekamp–Massey to get the shortest recurrence of $(s_k)$: with high probability it is the minimal polynomial $m(x)=x^d+\cdots +m_{1}x+m_0$ of $A$ (or of $A$ restricted to the Krylov space of $v$). Then: solve $Ax=b$ with $v=b$: if $m_0\ne 0$, $x=-\frac{1}{m_0}{\sum }_{i\ge 1}{m}_i{A}^{i-1}b$; determinant: precondition with a random diagonal $D$ so that the minimal and characteristic polynomials of $AD$ coincide, then $\det A=(-1{)}^n{m}_0/\det D$; rank similarly with preconditioners. For $n=10^6$ and $\mathrm{nnz}=10^7$ this is feasible where elimination is not. Uses: the linear algebra stage of the quadratic and number field sieves and of index-calculus discrete logs (millions of rows over ${\mathbb{F}}_2$, block variants), sparse determinants in computer algebra. Block Wiedemann (Coppersmith) uses blocks of $b$ vectors: matrix–matrix products with better cache behaviour that parallelise across machines.

vector<mint> berlekampMassey(const vector<mint>& s) {      // shortest c with s[i] = sum_j c[j] s[i-1-j]
    int n = sz(s), L = 0, m = 0; vector<mint> C(n + 1, 0), B(n + 1, 0), T; C[0] = B[0] = 1; mint b = 1;
    for (int i = 0; i < n; i++) {
        m++; mint d = s[i]; for (int j = 1; j <= L; j++) d += C[j] * s[i - j];
        if (d == mint(0)) continue;
        T = C; mint coef = d / b;
        for (int j = m; j <= n; j++) C[j] -= coef * B[j - m];
        if (2 * L > i) continue;
        L = i + 1 - L; B = T; b = d; m = 0;
    }
    vector<mint> c(L); for (int j = 0; j < L; j++) c[j] = -C[j + 1];
    return c;
}
 
struct SparseMat {                                          // n x n, list of (row, col, value)
    int n; vector<tuple<int, int, mint>> e;
    vector<mint> mul(const vector<mint>& v) const { vector<mint> r(n, 0); for (auto& [i, j, x] : e) r[i] += x * v[j]; return r; }
};
 
mt19937_64 wrng(12345);
vector<mint> wiedemannPoly(const SparseMat& A, const vector<mint>& u, vector<mint> v) {   // minimal polynomial of (u^T A^k v), monic, low degree first
    int n = A.n; vector<mint> s(2 * n);
    for (int k = 0; k < 2 * n; k++) { mint t = 0; for (int i = 0; i < n; i++) t += u[i] * v[i]; s[k] = t; v = A.mul(v); }
    auto c = berlekampMassey(s); int L = sz(c); vector<mint> f(L + 1); f[L] = 1;
    for (int j = 0; j < L; j++) f[L - 1 - j] = -c[j];
    return f;
}
 
vector<mint> sparseSolve(const SparseMat& A, const vector<mint>& b, int tries = 3) {   // nonsingular A; {} on failure
    int n = A.n;
    for (int t = 0; t < tries; t++) {
        vector<mint> u(n); for (auto& x : u) x = wrng() % 998244353;
        auto f = wiedemannPoly(A, u, b);
        if (sz(f) < 2 || f[0] == mint(0)) continue;
        vector<mint> x(n, 0), cur = b;                     // x = -(1/f0) sum_{i>=1} f_i A^{i-1} b
        for (int i = 1; i < sz(f); i++) { for (int j = 0; j < n; j++) x[j] += f[i] * cur[j]; cur = A.mul(cur); }
        mint inv = -f[0].inv(); for (auto& y : x) y *= inv;
        if (A.mul(x) == b) return x;
    }
    return {};
}
 
mint sparseDet(const SparseMat& A, int tries = 3) {         // Las Vegas for nonsingular A; returns 0 if every try fails (singular w.h.p.)
    int n = A.n; if (n == 0) return 1;
    for (int t = 0; t < tries; t++) {
        vector<mint> d(n); for (auto& x : d) x = wrng() % 998244352 + 1;
        SparseMat B = A; for (auto& [i, j, x] : B.e) x *= d[j];   // A D
        vector<mint> u(n), v(n); for (auto& x : u) x = wrng() % 998244353; for (auto& x : v) x = wrng() % 998244353;
        auto f = wiedemannPoly(B, u, v);
        if (sz(f) != n + 1) continue;
        mint det = n % 2 ? -f[0] : f[0], dd = 1; for (auto& x : d) dd *= x;
        return det / dd;
    }
    return 0;
}

Proof: Wiedemann

The sequence satisfies $m$. If $m(A)={\sum }_i{m}_i{A}^i=0$ then for every $k$, ${\sum }_i{m}_i{s}_{k+i}=u^T{A}^{k}m(A)v=0$, so $(s_k)$ satisfies the recurrence of the minimal polynomial; Berlekamp–Massey returns the shortest recurrence of the first $2n$ terms, which is the true minimal recurrence of the sequence because any two recurrences of length $\le n$ that agree on $2n$ terms agree forever. That minimal recurrence divides $m$; it is a proper divisor only if $u$ or $v$ is unluckily aligned with an invariant subspace, which by the Schwartz–Zippel lemma has probability $O(n/q)$ over a field of size $q$ (so a large prime makes one try almost always enough). Solve: with $v=b$, the recurrence polynomial $f$ (w.h.p. the minimal polynomial of $b$ under $A$) gives $f(A)b=0$, i.e. $f_{0}b=-A{\sum }_{i\ge 1}{f}_i{A}^{i-1}b$, and $f_0\ne 0$ for nonsingular $A$; the answer is verified, so the method is Las Vegas. Determinant: for nonsingular $A$ and a random diagonal $D$, the characteristic polynomial of $AD$ is squarefree with high probability (Kaltofen–Saunders / Chen et al. diagonal preconditioning), so a minimal polynomial of degree $n$ equals it and $\det (AD)=(-1{)}^n{f}_0$; divide by $\det D$. A degree-$n$ result is the characteristic polynomial, so it is always correct (for singular $A$ it has $f_0=\pm \det (AD)=0$ and the function returns 0); if no try reaches degree $n$ the function returns 0, which is right for singular $A$ and a low-probability failure otherwise. Berlekamp–Massey is proved in Linear Recurrences. Checked against dense elimination on random sparse matrices.

Structured matrices

template <class T> vector<T> polyStrip(vector<T> a) { while (!a.empty() && a.back() == T(0)) a.pop_back(); return a; }
template <class T> vector<T> polyMod(vector<T> a, vector<T> b) {
    a = polyStrip(a), b = polyStrip(b); assert(!b.empty());
    T inv = T(1) / b.back(); int m = sz(b) - 1;
    for (int n = sz(a) - 1; n >= m; n--) { T x = a[n] * inv; if (x == T(0)) continue; for (int j = 0; j <= m; j++) a[n - m + j] -= x * b[j]; }
    a.resize(min<int>(sz(a), m)); return polyStrip(a);
}
template <class T> T resultant(vector<T> a, vector<T> b) {
    a = polyStrip(a), b = polyStrip(b);
    if (a.empty() || b.empty()) return T(0);
    T res = 1;
    while (true) {
        if (sz(b) == 1) { for (int i = 0; i + 1 < sz(a); i++) res *= b[0]; return res; }   // Res(A, c) = c^deg A
        auto c = polyMod(a, b); if (c.empty()) return T(0);
        int da = sz(a) - 1, db = sz(b) - 1, dc = sz(c) - 1;
        if ((ll)da * db % 2) res = -res;
        for (int i = 0; i < da - dc; i++) res *= b.back();
        a = b, b = c;
    }
}
template <class T> T circulantDeterminant(const vector<T>& c) {     // M[i][j] = c[(j - i) mod n]
    int n = sz(c); if (n == 0) return T(1);
    auto C = polyStrip(c); if (C.empty()) return T(0);
    vector<T> A(n + 1, T(0)); A[0] = T(-1), A[n] = T(1);
    return resultant(A, C);
}
template <class T> T permutantDeterminant(const vector<T>& a, const vector<int>& p) {   // row 0 = a, M[i+1][j] = M[i][p[j]]
    int n = sz(a); if (n == 0) return T(1);
    vector<int> c(n, 0);
    for (int s = 1; s < n; s++) { c[s] = p[c[s - 1]]; if (c[s] == 0) return T(0); }   // more than one cycle: equal rows
    vector<T> b(n); for (int s = 0; s < n; s++) b[s] = a[c[s]];
    vector<char> vis(n, 0); int cycles = 0;
    for (int i = 0; i < n; i++) if (!vis[i]) { cycles++; for (int j = i; !vis[j]; j = c[j]) vis[j] = 1; }
    T d = circulantDeterminant(b);
    return ((n - cycles) + (n - 1) / 2) % 2 ? -d : d;
}
template <class T> T rankOneDeterminant(const vector<T>& a, const vector<T>& b, T x) {   // det(x I + a b^T)
    int n = sz(a); T s = 0, xp = 1;
    for (int i = 0; i < n; i++) s += a[i] * b[i];
    for (int i = 0; i + 1 < n; i++) xp *= x;
    return n == 0 ? T(1) : xp * x + xp * s;
}
template <class T> T vandermondeDeterminant(const vector<T>& x) {
    T d = 1; for (int i = 0; i < sz(x); i++) for (int j = i + 1; j < sz(x); j++) d *= x[j] - x[i]; return d;
}
template <class T> vector<vector<T>> vandermondeInverse(const vector<T>& x) {   // M[i][j] = x_i^j, x distinct
    int n = sz(x); vector<T> full(n + 1, T(0)); full[0] = 1;                  // prod (t - x_i)
    for (int i = 0; i < n; i++) for (int j = n; j >= 0; j--) { full[j] *= -x[i]; if (j) full[j] += full[j - 1]; }
    vector<vector<T>> inv(n, vector<T>(n));
    for (int i = 0; i < n; i++) {
        vector<T> up(n); T rem = full[n];                                       // full / (t - x_i), synthetic division
        for (int j = n - 1; j >= 0; j--) { up[j] = rem; rem = full[j] + rem * x[i]; }
        T down = 1; for (int j = 0; j < n; j++) if (j != i) down *= x[i] - x[j];
        down = T(1) / down;
        for (int j = 0; j < n; j++) inv[j][i] = up[j] * down;
    }
    return inv;
}
template <class T> vector<T> thomas(vector<T> lo, vector<T> di, vector<T> up, vector<T> rhs) {   // lo[i] x[i-1] + di[i] x[i] + up[i] x[i+1] = rhs[i]
    int n = sz(di);
    for (int i = 1; i < n; i++) { T w = lo[i] / di[i - 1]; di[i] -= w * up[i - 1]; rhs[i] -= w * rhs[i - 1]; }
    vector<T> x(n);
    for (int i = n - 1; i >= 0; i--) x[i] = (rhs[i] - (i + 1 < n ? up[i] * x[i + 1] : T(0))) / di[i];
    return x;
}
template <class T> bool shermanMorrison(vector<vector<T>>& B, const vector<T>& u, const vector<T>& v) {   // B = A^-1 becomes (A + u v^T)^-1
    int n = sz(B); vector<T> bu(n, T(0)), vb(n, T(0));
    for (int i = 0; i < n; i++) for (int j = 0; j < n; j++) bu[i] += B[i][j] * u[j], vb[j] += v[i] * B[i][j];
    T den = 1; for (int i = 0; i < n; i++) den += v[i] * bu[i];
    if (den == T(0)) return false;
    T inv = T(1) / den;
    for (int i = 0; i < n; i++) for (int j = 0; j < n; j++) B[i][j] -= bu[i] * vb[j] * inv;
    return true;
}

Proof: the structured formulas

Circulant. With $\omega$ a primitive $n$-th root of unity (in a splitting field), the vector $({\omega }^{jk}{)}_j$ is an eigenvector of $M$ with eigenvalue ${\sum }_d{c}_d{\omega }^{dk}=C({\omega }^k)$, and these $n$ vectors are independent (a Vandermonde matrix), so $\det M={\prod }_{k}C({\omega }^k)$. Since $x^n-1$ is monic with roots ${\omega }^k$, this product is $\mathrm{Res}(x^n-1,C)$, a polynomial identity valid over any field. Resultant algorithm: $\mathrm{Res}(A,B)=(-1{)}^{\mathrm{deg}A\mathrm{deg}B}\mathrm{Res}(B,A)$ and, for $A=QB+R$, $\mathrm{Res}(B,A)=\mathrm{lc}(B{)}^{\mathrm{deg}A-\mathrm{deg}R}\mathrm{Res}(B,R)$ (both follow from $\mathrm{Res}(B,A)=\mathrm{lc}(B{)}^{\mathrm{deg}A}{\prod }_{B(\beta )=0}A(\beta )$ and $A(\beta )=R(\beta )$); a constant $c$ gives $c^{\mathrm{deg}A}$, and a zero remainder means a common factor, resultant 0. Permutant: rows are $a\circ p^i$; if $p$ has a cycle not containing 0, the orbit of 0 returns to 0 before step $n$ and two rows coincide. Otherwise, with $c_s=p^s(0)$, permuting columns by $c$ and reversing the row order turns $M$ into the circulant of $b_s=a_{c_s}$; the sign is that of the column permutation ($(-1{)}^{n-\#\text{cycles}}$) times the reversal ($\lfloor \frac{n-1}{2}\rfloor$ swaps).

Rank one. $a{b}^T$ has rank $\le 1$, with eigenvalue $b^{T}a$ (eigenvector $a$) and $0$ on the $(n-1)$-dimensional kernel of $b^T$, so $\det (xI+a{b}^T)=x^{n-1}(x+b^{T}a)$. Vandermonde: the determinant is a polynomial of degree $\left(\genfrac{}{}{0ex}{}{n}{2}\right)$ vanishing whenever $x_i=x_j$, hence divisible by ${\prod }_{i<j}(x_j-x_i)$, with matching leading coefficient. The inverse: ${\sum }_j{M}_{kj}(M^{-1}{)}_{ji}={\sum }_j{x}_k^j{\ell }_{i,j}={\ell }_i(x_k)=[k=i]$ for the Lagrange basis polynomial ${\ell }_i={\prod }_{j\ne i}\frac{t-x_j}{x_i-x_j}$, whose coefficients the code computes by synthetic division of $\prod (t-x_j)$ by $t-x_i$. Thomas is Gaussian elimination without pivoting on a tridiagonal matrix (each step touches one row), stable for diagonally dominant systems; lo[0] and up[n-1] are ignored. Sherman–Morrison: multiply $(A+u{v}^T)(B-\frac{Bu{v}^{T}B}{1+v^{T}Bu})$ out with $AB=I$: $I+u{v}^{T}B-\frac{u{v}^{T}B+u(v^{T}Bu)v^{T}B}{1+v^{T}Bu}=I$; if $1+v^{T}Bu=0$ then $(A+u{v}^T)Bu=u(1+v^{T}Bu)=0$ with $Bu\ne 0$, so the update is singular. All verified against dense determinants and inverses.

Lanczos and Krylov methods

For a symmetric $A$ and a unit start vector $v_1$, the Krylov space ${\mathcal{K}}_k=\mathrm{span}\{v_1,A{v}_1,\dots ,A^{k-1}{v}_1\}$ gets an orthonormal basis from a three-term recurrence
${\beta }_{j+1}{v}_{j+1}=A{v}_j-{\alpha }_j{v}_j-{\beta }_j{v}_{j-1},{\alpha }_j=v_j^{T}A{v}_j,{\beta }_{j+1}=\Vert A{v}_j-{\alpha }_j{v}_j-{\beta }_j{v}_{j-1}\Vert ,$
and in that basis $A$ becomes the tridiagonal $T_k$ (diagonal $\alpha$, off-diagonal $\beta$), whose extreme eigenvalues approximate those of $A$ long before $k=n$ (often a few dozen steps for millions of rows). Cost $O(k\cdot \mathrm{nnz}+k^2)$, memory $O(n)$ for eigenvalues only.

pair<vector<double>, vector<double>> lanczos(int n, const function<vector<double>(const vector<double>&)>& mulA, int k, mt19937& rng) {
    vector<vector<double>> V; vector<double> al, be;        // full reorthogonalisation: stable, O(k^2 n)
    vector<double> v(n); for (auto& x : v) x = uniform_real_distribution<double>(-1, 1)(rng);
    double nr = 0; for (double x : v) nr += x * x; nr = sqrt(nr); for (auto& x : v) x /= nr;
    for (int j = 0; j < k; j++) {
        V.push_back(v); vector<double> w = mulA(v);
        double a = 0; for (int i = 0; i < n; i++) a += w[i] * v[i]; al.push_back(a);
        for (auto& q : V) { double d = 0; for (int i = 0; i < n; i++) d += w[i] * q[i]; for (int i = 0; i < n; i++) w[i] -= d * q[i]; }
        double b = 0; for (double x : w) b += x * x; b = sqrt(b);
        if (j + 1 == k || b < 1e-10) break;
        be.push_back(b); for (int i = 0; i < n; i++) v[i] = w[i] / b;
    }
    return {al, be};                                        // T: diagonal al, off-diagonal be
}
 
double tridiagLargestEigenvalue(const vector<double>& al, const vector<double>& be) {   // Sturm bisection
    int k = sz(al); double lo = 1e18, hi = -1e18;
    for (int i = 0; i < k; i++) { double r = (i ? fabs(be[i - 1]) : 0) + (i < sz(be) ? fabs(be[i]) : 0); lo = min(lo, al[i] - r); hi = max(hi, al[i] + r); }
    auto below = [&](double x) { int cnt = 0; double d = 1; for (int i = 0; i < k; i++) { d = al[i] - x - (i ? be[i - 1] * be[i - 1] / d : 0); if (d == 0) d = -1e-300; cnt += d < 0; } return cnt; };
    for (int it = 0; it < 200; it++) { double mid = (lo + hi) / 2; if (below(mid) == k) hi = mid; else lo = mid; }
    return hi;
}

	Lanczos	Wiedemann
needs symmetry	yes (use $A^{T}A$ otherwise)	no
iterations	$\approx n$	$\approx 2n$
failure mode	self-orthogonal vectors over ${\mathbb{F}}_2$ (handled by blocking)	unlucky random vectors (retry)
parallelism	good (block Lanczos, Montgomery)	excellent (block Wiedemann)
memory	$O(n)$ vectors	$O(n)$ scalars

Krylov method	For
Lanczos	symmetric eigenproblems
conjugate gradient	symmetric positive definite systems (Lanczos rearranged)
Arnoldi	non-symmetric eigenproblems, full orthogonalisation $O(k^{2}n)$
GMRES	non-symmetric systems, built on Arnoldi
power iteration	the dominant eigenvalue only; Lanczos is strictly better

In floating point the plain recurrence loses orthogonality once an eigenvalue converges and produces “ghost” copies; the code reorthogonalises fully (fine for small $k$), real implementations use selective or partial reorthogonalisation; over a finite field the issue disappears. Uses: spectral graph theory (second Laplacian eigenvalue, clustering, expansion), PageRank-type problems, block Lanczos over ${\mathbb{F}}_2$ for factorisation, ground states of sparse Hamiltonians.

Proof: the three-term recurrence and Sturm counting

By construction $v_{j+1}\perp v_j,v_{j-1}$. For $i<j-1$: $v_i^{T}A{v}_j=(A{v}_i{)}^T{v}_j$ (symmetry), and $A{v}_i\in \mathrm{span}(v_1..v_{i+1})$ with $i+1<j$, so it is orthogonal to $v_j$; hence $A{v}_j$ has no component along older vectors and the recurrence produces an orthonormal basis $V$ with $V^{T}AV=T$ tridiagonal. If $k$ reaches the dimension of the Krylov space, $\mathcal{K}$ is $A$-invariant and the eigenvalues of $T$ are eigenvalues of $A$. Sturm count: for symmetric tridiagonal $T-xI=LD{L}^T$, the pivots satisfy $d_i={\alpha }_i-x-{\beta }_{i-1}^2/d_{i-1}$, and by Sylvester’s law of inertia the number of negative pivots equals the number of eigenvalues below $x$; bisection on that count isolates the largest eigenvalue (Gershgorin bounds start the interval). Checked on $A=QD{Q}^T$ with random orthogonal $Q$ and known $D$.

Related notes

• Linear Recurrences (Berlekamp–Massey)
• Polynomial GCD and Resultant
• Gaussian Elimination