Fraction-free elimination and echelon forms

Gaussian elimination over \mathbb{Q} applied to an integer matrix creates fractions whose numerators and denominators grow quickly. Fraction-free elimination avoids them: it multiplies rows by pivots instead of dividing by them, and divides out the previous pivot exactly, so that every intermediate entry is a minor of the input. fmpz_mat_fflu computes such a decomposition, from which determinants, ranks and solutions follow; fmpz_mat_rref computes the reduced row echelon form, either by fraction-free elimination or by a multimodular method with an exact certificate.

Fraction-free LU decomposition#

fmpz_mat_fflu(B, den, perm, A, rank_check) performs Bareiss' algorithm [Bar1968] on the m \times n matrix A and returns its rank r. Pivots are chosen with fmpz_mat_find_pivot_any: in each column, the first row at or below the current pivot row with a nonzero entry; rows are swapped to bring it into place, and the swaps are applied to perm if it is not NULL. When the k-th pivot p_k has been chosen, every entry of the remaining submatrix is updated as

with p_{-1} = 1, where c_k is the pivot column; the division is exact. After the elimination:

• the first r rows of B, from their pivot columns on, form the fraction-free echelon form: the pivot in row k is p_k, the determinant of the leading (k+1) \times (k+1) submatrix of the row-permuted A restricted to the first k + 1 pivot columns;
• below each pivot, B keeps the entries of the pivot column as they were when the pivot was chosen (the multipliers of the fraction-free LU decomposition [NakTurWil1997]); the other entries of the rows below the rank are zero;
• den is the last pivot p_{r-1}, or 1 if r = 0.

For a nonsingular square A, den is \det(PA) = \pm \det(A), where P is the row permutation. With rank_check set, the function stops at the first column without pivot and returns 0 with den = 0, which is what the determinant and solving functions need. B may be A.

The entries of B are minors of A, so they have at most about n(b + \frac12 \log_2 n) bits for b-bit entries, and the cost is O(mnr) multiplications and exact divisions of such integers. For small entries (at most 62 bits on a 64-bit machine) the updates use two-word arithmetic and division by a precomputed inverse of the previous pivot until an entry outgrows a word.

Reduced row echelon form#

fmpz_mat_rref(B, den, A) sets (B, den) to the reduced row echelon form of A: B/\mathit{den} is in reduced row echelon form, so every pivot entry of B equals den, the other entries of the pivot columns are zero, and the rows below the rank are zero. It returns the rank. The denominator is not normalised: it is a nonzero r \times r minor of A for fmpz_mat_rref_fflu, and 1 or a divisor of such a minor for fmpz_mat_rref_mul, so the two can return different but equivalent results. Divide by the gcd of den and the entries for a canonical form, or use fmpq_mat_rref.

For an m \times n matrix, fmpz_mat_rref calls fmpz_mat_rref_fflu if m \le 3, if n \le 2, if m \le 20 and n > m, or if 20 < m \le 100 and n > m + \frac{m - 20}{80} m; otherwise it calls fmpz_mat_rref_mul. Wide matrices thus go to fraction-free elimination, whose cost is dominated by the elimination of the first m columns, and tall or square ones to the multimodular method.

Fraction-free back substitution#

fmpz_mat_rref_fflu computes the fraction-free echelon form with fmpz_mat_fflu. Writing its first r rows, up to a permutation of the columns, as (U \mid V) with U upper triangular, the reduced form is (I \mid U^{-1} V), and \mathit{den} \cdot U^{-1} V is integral since \mathit{den} = \det U up to sign. Its columns are computed by back substitution with exact divisions by the pivots, scaled to the common denominator den; only the non-pivot columns need to be updated. This is the fraction-free Gauss–Jordan elimination of [NakTurWil1997] restricted to the part V.

Multimodular echelon form#

fmpz_mat_rref_mul follows [Stein2007]:

1. Reduce A modulo a prime p of 17 bits and compute its rank r, pivot rows and pivot columns modulo p with nmod_mat_lu_with_pivots. If r = n, the reduced form is the identity matrix on top of zero rows, and \mathit{den} = 1.
2. Let B_0 be the r \times r submatrix of A formed by the pivot rows and columns, which is nonsingular modulo p and hence over \mathbb{Q}, and C the pivot rows restricted to the other columns. Solve B_0 E_2 = \mathit{den} \cdot C with fmpz_mat_solve; the candidate is \mathit{den} on the pivot positions and E_2 in the other columns.
3. Certify: check that the candidate is in reduced row echelon form, and that every non-pivot row of A is annihilated by the nullspace basis defined by the candidate (one matrix product). Both checks are exact, so if they pass, the rows of A lie in the span of the pivot rows, the rank is r and the candidate is the reduced row echelon form (fmpz_mat_rref_upper_certify_lu_mod_p).
4. If a check fails, the rank modulo p was below the rank of A; try the next prime.

The work is one solve with n - r right-hand sides, done by the output-sensitive modular solvers, and one product, instead of an elimination with growing entries.

Example#

A fraction-free LU decomposition of a 4 \times 4 matrix, whose pivots are the leading principal minors 2, 5, 18, 82, and the reduced row echelon form of a 3 \times 5 matrix of rank 2 computed in two ways:

#include <flint/flint.h>
#include <flint/fmpz.h>
#include <flint/fmpz_mat.h>

int main(void)
{
    fmpz_mat_t A, B;
    fmpz_t den;
    slong i, j, rank, perm[4] = { 0, 1, 2, 3 };

    fmpz_init(den);
    fmpz_mat_init(A, 4, 4);
    fmpz_mat_init(B, 4, 4);

    /* tridiagonal with i + 2 on the diagonal, plus A[0][3] = 3 */
    for (i = 0; i < 4; i++)
        for (j = 0; j < 4; j++)
            fmpz_set_si(fmpz_mat_entry(A, i, j),
                (i == j) ? i + 2 : (j == i + 1 || i == j + 1));
    fmpz_set_si(fmpz_mat_entry(A, 0, 3), 3);

    rank = fmpz_mat_fflu(B, den, perm, A, 0);
    flint_printf("rank %wd, den %{fmpz}\n", rank, den);
    fmpz_mat_print_pretty(B);
    flint_printf("\n");

    fmpz_mat_clear(A);
    fmpz_mat_clear(B);

    /* row 2 = row 0 + 2 row 1 */
    fmpz_mat_init(A, 3, 5);
    fmpz_mat_init(B, 3, 5);
    for (j = 0; j < 5; j++)
    {
        fmpz_set_si(fmpz_mat_entry(A, 0, j), 2 * j + 2);
        fmpz_set_si(fmpz_mat_entry(A, 1, j), j * j - 3);
        fmpz_set_si(fmpz_mat_entry(A, 2, j), 2 * j + 2 + 2 * (j * j - 3));
    }

    rank = fmpz_mat_rref_fflu(B, den, A);
    flint_printf("rref_fflu: rank %wd, den %{fmpz}\n", rank, den);
    fmpz_mat_print_pretty(B);
    flint_printf("\n");

    rank = fmpz_mat_rref_mul(B, den, A);
    flint_printf("rref_mul: rank %wd, den %{fmpz}\n", rank, den);
    fmpz_mat_print_pretty(B);
    flint_printf("\n");

    fmpz_mat_clear(A);
    fmpz_mat_clear(B);
    fmpz_clear(den);
    return 0;
}

rank 4, den 82
[[2 1 0 3]
[1 5 2 -3]
[0 2 18 8]
[0 0 5 82]
]
rref_fflu: rank 2, den 8
[[8 0 -16 -40 -72]
[0 8 20 36 56]
[0 0 0 0 0]
]
rref_mul: rank 2, den 8
[[8 0 -16 -40 -72]
[0 8 20 36 56]
[0 0 0 0 0]
]

Both echelon forms have the denominator 8, while the reduced fraction form has denominator 2 since all entries are divisible by 4. In the decomposition, the entries 1, 2 and 5 below the diagonal are the multipliers.