Engineering/Mathematics/Algebra and algorithms

Linearly Generated Sequences and Sparse Linear Systems

A sequence satisfying a linear recurrence is completely described by one polynomial, and that polynomial can be recovered from twice its degree in terms. Applied to the Krylov sequence of a sparse matrix, this turns an intractable dense elimination into a linear number of matrix–vector products.

Advanced theory
Algebra and algorithms
Sieve back-end
≈15 min read
Enables record factorizations

2n termsRecovery thresholdA sequence with a minimal polynomial of degree ≤ n is determined by its first 2n terms.

O(n·w)Wiedemann costn matrix–vector products at w non-zeros each, versus O(n³) for dense elimination.

gcdOne algorithmThe extended Euclidean algorithm computes minimal polynomials, decodes Reed–Solomon codes and performs rational reconstruction — the same procedure each time.

No fill-inMemory advantageIterative methods never modify the matrix, so sparsity is preserved throughout.

Executive summary

A sequence (a₀, a₁, …) over a field is linearly generated if some non-zero polynomial g(X) = ∑ c_iXⁱ annihilates it, meaning ∑_i c_i a_j+i = 0 for all j. The set of such polynomials is an ideal, so it has a unique monic generator: the minimal polynomial of the sequence.

Two facts make this useful. The minimal polynomial can be computed from 2n terms by the extended Euclidean algorithm — equivalently by Berlekamp–Massey — in O(n²) field operations. And for a matrix A and vectors u, v, the scalar sequence a_i = u^TAⁱv is linearly generated by the minimal polynomial of A. Wiedemann's algorithm combines the two to solve sparse systems using only matrix–vector products.

ObjectLinearly generated sequence

Annihilated by a polynomial; the annihilators form an ideal in F[X].

InvariantMinimal polynomial

The monic generator; degree n means 2n terms suffice to determine it.

AlgorithmExtended Euclid / Berlekamp–Massey

Recovers the minimal polynomial in O(n²) operations.

ApplicationWiedemann

Solve Ax = b or find ker A using n matrix–vector products and no fill-in.

Contents

Linearly generated sequences

Definition D1

Generator and minimal polynomial

For an infinite sequence (a_i) over a field F, a polynomial g = ∑_i=0^k c_iXⁱ generates the sequence if ∑_i=0^k c_ia_j+i = 0 for every j ≥ 0. The generators together with 0 form an ideal of F[X]; if it is non-zero its monic generator is the minimal polynomial of the sequence, and it divides every other generator.

Degree meaning. A minimal polynomial of degree k corresponds to the shortest linear feedback shift register producing the sequence, and its length is exactly the linear complexity of the sequence.
Determination threshold. If the minimal polynomial has degree at most n, the first 2n terms determine it uniquely. This is the same information-theoretic count that appears in Reed–Solomon decoding and in rational reconstruction — because they are the same computation.
Matrix case. For a square matrix A, the minimal polynomial of the sequence a_i = u^TAⁱv divides the minimal polynomial of A, which in turn divides the characteristic polynomial. For random u, v the divisor is usually the full minimal polynomial of A.

∑_{i ≥ 0} a_iXⁱ = f(X)/g^rev(X)A sequence is linearly generated exactly when its generating series is a rational function — the link that makes rational function reconstruction the natural algorithm.

Contents

Computing the minimal polynomial

Given 2n terms, the minimal polynomial is recovered by running the extended Euclidean algorithm on X²ⁿ and the polynomial formed from the terms, stopping when the remainder degree drops below n.

Minimal polynomial from 2n terms

input: a₀ … a_{2n−1} ∈ F
A(X) ← ∑_{i<2n} a_i X^{2n−1−i}
run extended Euclid on (X^{2n}, A), keeping the coefficient sequence
stop at the first remainder r with deg r < n
the corresponding multiplier t gives the reversed minimal polynomial
normalise t to monic and reverse its coefficients

O(n²) field operations classically, O(n log² n) with fast polynomial arithmetic. Berlekamp–Massey performs the same computation in an incremental form with the same asymptotic cost.

Three views of the same procedure
Problem	Input	Output	Stopping rule
Minimal polynomial	2n sequence terms	shortest recurrence	remainder degree < n
Rational function reconstruction	z and modulus f	r/t with r ≡ zt	degree bounds on r and t
Reed–Solomon decoding	received word syndromes	error locator polynomial	half the redundancy consumed
Rational reconstruction over ℤ	residue and modulus	rational number	size bounds on numerator and denominator

Recognising these as one algorithm is a genuine simplification: a single well-tested extended-gcd routine with configurable stopping conditions serves all four.

Contents

Wiedemann's algorithm for sparse systems

The matrices produced by sieve algorithms have millions of rows but only tens of non-zero entries per row. Dense elimination is impossible; iterative methods that touch the matrix only through multiplication are the practical answer.

Form the Krylov sequence
Choose random vectors u, v and compute a_i = u^TAⁱv for i = 0, …, 2n−1 using 2n matrix–vector products.
Recover the minimal polynomial
Run the extended Euclidean procedure on the 2n scalars to obtain g, the minimal polynomial of the sequence — with high probability the minimal polynomial of A itself.
Use g to solve the problem
If g(0) ≠ 0, the identity g(A) = 0 expresses A⁻¹ as a polynomial in A, giving a solution of Ax = b. If g(0) = 0, factor out X^m to obtain a kernel vector.
Verify and repeat if necessary
Check the candidate directly with one more matrix–vector product. Failure occurs only for unlucky u, v; a fresh random choice succeeds with high probability.

Cost comparison on an n × n matrix with w non-zeros per row
Method	Time	Space	Suitability
Dense Gaussian elimination	O(n³)	O(n²)	Small dense systems only
Structured Gaussian elimination	problem-dependent	grows with fill-in	Pre-processing stage to shrink the matrix
Wiedemann	O(n²w) total	O(nw)	Huge sparse systems; simple to implement
Block Wiedemann / block Lanczos	same order, K-way parallel	O(nw)	Production factoring records; distributed across machines

The Wiedemann time is n matrix–vector products of cost O(nw) each, plus O(n²) for the sequence step. Over F₂ the vector operations are word-parallel XORs, which is why the block variants pack 64 or 128 right-hand sides into machine words.

Why the block variants exist

Blocking replaces the scalar sequence by a matrix sequence, using K random vectors at once. The number of iterations drops by a factor of K, the work per iteration is unchanged in word terms because K vectors fit in one machine word per coordinate, and the iterations parallelise cleanly. This is what makes the linear algebra stage of a large factoring effort feasible on a cluster.

Contents

Other uses of the same machinery

Minimal polynomial of a field element

Coordinates of the powers of α in F[X]/(f) form a linearly generated sequence; its minimal polynomial is the minimal polynomial of α.

Stream cipher analysis

Berlekamp–Massey recovers the shortest LFSR from a keystream segment, which is why linear complexity is a basic design criterion for keystream generators.

Reed–Solomon decoding

The syndrome sequence is linearly generated by the error locator polynomial; recovering it is exactly this computation.

Sequence prediction

Any sequence satisfying a fixed-order linear recurrence — including combinatorial ones like Fibonacci-type sequences — is identified from a finite prefix.

Order of a matrix

The minimal polynomial gives the eigenstructure over an extension field and hence the multiplicative order of an invertible matrix.

Discrete log linear algebra

Index-calculus systems over ℤ_q are solved by the same iterative approach, with modular rather than binary arithmetic.

Contents

Quick reference and FAQ

Facts
Fact	Statement
Ideal structure	Generators of a sequence form an ideal of F[X]
Determination	Degree ≤ n minimal polynomial is fixed by 2n terms
Rationality	Linearly generated ⟺ generating series is rational
Krylov sequence	u^TAⁱv is generated by the minimal polynomial of A
Cayley–Hamilton	Minimal polynomial divides the characteristic polynomial
Wiedemann iterations	2n matrix–vector products, plus O(n²) for the sequence step
Kernel extraction	g(0) = 0 signals a singular matrix and yields a kernel vector

Is Berlekamp–Massey different from the extended Euclidean algorithm?

They compute the same thing and have the same complexity. Berlekamp–Massey is incremental — it processes one term at a time and is convenient in streaming or hardware settings — while the Euclidean formulation reuses an existing gcd routine and generalises directly to rational function reconstruction.

What is the failure probability of Wiedemann?

For random u, v over a field of size q, the recovered polynomial is the true minimal polynomial with probability at least about 1 − 2n/q. Over F₂ that bound is useless, which is precisely why block variants with many random vectors are used there; failure is detected cheaply by verification and handled by retry.

When is structured Gaussian elimination still worth doing?

Almost always as a pre-processing step. Sieve matrices contain many columns with a single non-zero entry — relations involving a large prime that appears once — and eliminating those shrinks the matrix substantially at negligible fill-in cost before the iterative phase begins.

Why not use a numerical iterative solver such as conjugate gradient?

Those methods rely on an inner product being positive definite, which is meaningless over a finite field. Lanczos can be adapted, but requires care because self-orthogonal vectors exist over F₂; the block Lanczos variant handles this explicitly.

Contents

References and further reading

V. Shoup, A Computational Introduction to Number Theory and Algebra, Cambridge University Press, 2005 — Chapter 19.
D. H. Wiedemann, 'Solving sparse linear equations over finite fields', IEEE Trans. Inform. Theory 32 (1986) 54–62.
P. L. Montgomery, 'A block Lanczos algorithm for finding dependencies over GF(2)', EUROCRYPT '95, LNCS 921, 106–120.
J. L. Massey, 'Shift-register synthesis and BCH decoding', IEEE Trans. Inform. Theory 15 (1969) 122–127.

KEVOS® Knowledge LibraryEngineering → MathematicsTaxonomy ID: ENG-MATHPage ID: sparse-linear-systems-and-wiedemannReview cycle: annual

Linearly Generated Sequences and Sparse Linear Systems

Executive summary

Linearly generated sequences

Generator and minimal polynomial

Computing the minimal polynomial

Wiedemann's algorithm for sparse systems

Form the Krylov sequence

Recover the minimal polynomial

Use g to solve the problem

Verify and repeat if necessary

Other uses of the same machinery

Minimal polynomial of a field element

Stream cipher analysis

Reed–Solomon decoding

Sequence prediction

Order of a matrix

Discrete log linear algebra

Quick reference and FAQ

References and further reading

Continue learning

KEVOS Intelligence

Executive summary

Linearly generated sequences

Generator and minimal polynomial

Computing the minimal polynomial

Wiedemann's algorithm for sparse systems

Form the Krylov sequence

Recover the minimal polynomial

Use g to solve the problem

Verify and repeat if necessary

Other uses of the same machinery

Minimal polynomial of a field element

Stream cipher analysis

Reed–Solomon decoding

Sequence prediction

Order of a matrix

Discrete log linear algebra

Quick reference and FAQ

Related pages in this library

References and further reading

Continue learning