Fraction-free elimination

Let R be an integral domain, such as \mathbb{Z}, \mathbb{Z}[x] or \mathbb{F}_p[x], and A an m \times n matrix over R. Gaussian elimination over the fraction field creates quotients whose numerators and denominators grow; removing common factors needs gcd computations, which are expensive or unavailable. Fraction-free elimination keeps all entries in R, uses only exact divisions, and bounds every intermediate entry by a minor of A [Bar1968; NakTurWil1997]. It is the method of choice for small matrices over \mathbb{Z} and for matrices over polynomial rings, where modular methods are less developed.

Generic version#

gr_mat_fflu implements the algorithm below over any ring with exact division, and gr_mat_det_fflu, gr_mat_rank_fflu, gr_mat_nonsingular_solve_fflu, gr_mat_nonsingular_solve_den_fflu and gr_mat_rref_den_fflu build on it. The specialised versions are separate implementations with the same structure:

• fmpz_mat_fflu over \mathbb{Z} chooses the first nonzero pivot, and while the entries fit in a word computes the updates in two-word arithmetic with a precomputed inverse of the previous pivot;
• fmpz_poly_mat_fflu and nmod_poly_mat_fflu over \mathbb{Z}[x] and (\mathbb{Z}/n\mathbb{Z})[x] choose, among the nonzero entries of the column, one of least degree and, over \mathbb{Z}[x], with the smallest coefficients (fmpz_poly_mat_find_pivot_partial), which keeps the entries small;
• ca_mat_fflu over exact real and complex numbers chooses pivots with zero tests that may fail (linear algebra over exact numbers).

Bareiss' algorithm#

The divisions in step 2.3 are exact. This follows from Sylvester's identity.

With row exchanges, the same holds for the row-permuted matrix, and with skipped columns for the matrix restricted to the pivot columns and the column j. Consequently every intermediate entry is a minor of A: over \mathbb{Z} it is bounded by Hadamard's inequality (Bounds), with about n(b + \frac12 \log_2 n) bits for b-bit entries, and over K[x] its degree is at most n times the largest degree. Gaussian elimination over \mathbb{Q} without gcds would instead produce entries whose sizes grow exponentially with k.

For a nonsingular n \times n matrix, the last pivot is p_{n-1} = \det(PA) = \operatorname{sgn}(P) \det A, so fmpz_mat_det_bareiss returns den with the sign of the permutation. Each of the O(mnr) updates costs two multiplications and one exact division of minors, so over \mathbb{Z} the cost is O(n^3 M(n(b + \log n))) word operations for an n \times n matrix, that is O(n^5 b^2) with classical integer multiplication. This is competitive with the modular methods of Multimodular and p-adic linear algebra over the integers only for small n (FLINT's crossovers are n \approx 15 to 25).

The fraction-free LU decomposition#

FLINT stores the result in place, like an LU decomposition: the upper part of B, from the pivot of each row on, holds the fraction-free echelon form U, and below each pivot p_k the entries of column c_k are kept as they were when p_k was chosen (step 2.3 does not change column c). Let L be the lower triangular matrix formed by these columns, with the pivots p_k on its diagonal.

This is the fraction-free LU decomposition of [NakTurWil1997; Jef2010]. For the 4 \times 4 example of Example, with pivots 2, 5, 18, 82, one checks A = L \operatorname{diag}(2, 10, 90, 1476)^{-1} U.

Solving#

To solve AX = B for nonsingular A, fmpz_mat_solve_fflu_precomp applies to PB the same fraction-free updates as step 2.3 (forward substitution: the column j of the right-hand side is updated like a column of A), which yields the last column of the elimination of the augmented matrix (A \mid B), and then performs back substitution

where y is the transformed right-hand side. The result is X = p_{n-1} A^{-1} B, which is integral by Cramer's rule since p_{n-1} = \pm \det A, so all divisions are exact. fmpz_mat_solve_fflu then fixes the sign, returning \mathit{den} = \det A. For a rank-deficient A the same procedure, with the non-pivot unknowns set to zero, either produces a candidate or meets an inexact division and reports failure; fmpz_mat_can_solve_fflu then checks the remaining equations exactly.

Reduced echelon form#

Let the first r rows of the echelon form be (U \mid V) up to a permutation of the columns, where U is the r \times r upper triangular matrix of pivot columns, with the pivots p_0, \ldots, p_{r-1} on its diagonal. The reduced echelon form over the fraction field is (I \mid U^{-1} V) in the original column order. Since (U \mid V) is obtained from the pivot rows (A_P \mid A_N) of A by invertible row operations, U^{-1} V = A_P^{-1} A_N, and p_{r-1} U^{-1} V = \pm \det(A_P) A_P^{-1} A_N is integral by Cramer's rule. fmpz_mat_rref_fflu and fmpz_poly_mat_rref compute it by back substitution with exact divisions by the pivots, one non-pivot column at a time, and set the pivot columns to p_{r-1} times unit vectors. This is fraction-free Gauss–Jordan elimination [NakTurWil1997] restricted to the non-pivot columns, which are all that change.

Use in FLINT#