Computing an LU decomposition

An LU decomposition with row pivoting writes an $m \times n$ matrix $A$ over a field $\mathbb{Z}/p\mathbb{Z}$ as a row permutation of a product $LU$, with $L$ lower triangular and $U$ in row echelon form. It is the basis of nmod_mat_solve, nmod_mat_inv, nmod_mat_det, nmod_mat_rank, nmod_mat_rref and nmod_mat_nullspace, and it is useful directly when several systems with the same matrix must be solved. The functions overwrite $A$ with both factors and return the rank.

Output#

Let $r$ be the returned rank. On return:

• P[i] is the index of the row of the original $A$ that ended up in row $i$. P must have room for $m$ entries and need not be initialised.
• $U$ is the $r \times n$ matrix formed by the first $r$ rows of the output, keeping the entries in columns $j \ge i$ of row $i$ and replacing those in columns $j < i$ by zero. It is in row echelon form: the pivot (first nonzero entry) of row $i$ is in column $c_i$, with $c_0 < c_1 < \cdots < c_{r-1}$.
• $L$ is the $m \times r$ unit lower triangular matrix whose entry $(i, j)$ for $j < \min(i, r)$ is entry $(i, j)$ of the output; its diagonal entries are 1 and are not stored. The rows $i \ge r$ of the output contain $L$ in their first $r$ columns and zeros elsewhere.

Then row $i$ of $LU$ equals row P[i] of the original $A$. In matrix terms, $A = \Pi L U$ where $\Pi$ is the permutation matrix with a 1 in row P[i] and column $i$. nmod_mat_permute_rows(A, P, NULL) applies the same permutation to another copy of $A$, as in the example below.

The pivots are chosen column by column: the pivot of column $c$ is the first row, among those not yet used, with a nonzero entry in column $c$ after elimination. The output therefore does not depend on the algorithm, and $c_0, \ldots, c_{r-1}$ is the column rank profile of $A$, the lexicographically smallest set of $r$ linearly independent columns. nmod_mat_lu_with_pivots returns these indices, followed by the indices of the other columns.

Checking the rank#

With a nonzero rank_check, the decomposition stops as soon as a column without a pivot is found, and the function returns 0, leaving $A$ in an undefined state. This is what nmod_mat_solve and nmod_mat_inv use to reject singular matrices early. If $m \ge n$, the return value is then 0 exactly when the rank is less than $n$. If $m < n$ the result depends on the algorithm: for example, nmod_mat_lu_recursive returns 0 for the $2 \times 6$ matrix $(I_2 \mid 0)$, which has full rank 2. Use rank_check only with $m \ge n$.

Moduli#

The algorithms divide by pivots, so the modulus should be prime. For a composite modulus the functions still return a decomposition if every pivot they meet happens to be invertible, but they raise a FLINT_IMPINV exception when they meet a nonzero pivot that is not, and the returned "rank" has no meaning. Use nmod_mat_howell_form or nmod_mat_det for composite moduli.

Algorithms#

nmod_mat_lu_classical is Gaussian elimination, with one reduction for each update of an entry. nmod_mat_lu_classical_delayed performs the same operations but accumulates the updates of each entry of the active submatrix unreduced in one, two or three words, as in a dot product, and reduces an entry only when its column becomes the pivot column or its row the pivot row. nmod_mat_lu_recursive splits the columns into halves of $\lfloor n/2 \rfloor$ and $\lceil n/2 \rceil$, decomposes the left half recursively, updates the right half with a triangular solve (nmod_mat_solve_tril) and a matrix product (nmod_mat_submul), and decomposes the remaining block recursively. Its cost is dominated by matrix multiplication, $O(m n r^{\omega - 2})$ operations with the exponent $\omega$ of the multiplication algorithm used.

nmod_mat_lu uses the classical algorithm when $\min(m, n) \le 3$, or $\le 7$ for $p \ge 2^{28}$. Otherwise it uses the delayed algorithm when $\min(m, n)$ is below a crossover that depends on the bit size $b$ of $p$, and the recursive algorithm above it. The recursive algorithm calls nmod_mat_lu for the halves, so the crossover also sets its base case. The values were tuned on an AMD Zen 3 processor:

The details, with proofs and complexity bounds, are in Gaussian elimination and LU decomposition.

Example#

The decomposition of a $4 \times 5$ matrix of rank 3 modulo 17, and a check that $LU$ is $A$ with its rows permuted by P:

#include <flint/flint.h>
#include <flint/nmod_mat.h>

int main(void)
{
    /* a 4 x 5 matrix of rank 3 modulo 17: row 3 = row 0 + row 1 */
    ulong a[4][5] = { { 0, 0, 3, 1, 4 },
                      { 2, 4, 1, 0, 5 },
                      { 1, 2, 7, 7, 1 },
                      { 2, 4, 4, 1, 9 } };
    nmod_mat_t A, LU, L, U, T;
    slong P[4], piv[5], i, j, r;

    nmod_mat_init(A, 4, 5, 17);
    for (i = 0; i < 4; i++)
        for (j = 0; j < 5; j++)
            nmod_mat_entry(A, i, j) = a[i][j];

    nmod_mat_init_set(LU, A);
    r = nmod_mat_lu_with_pivots(P, piv, LU);

    flint_printf("rank %wd, P = [%wd %wd %wd %wd]\n", r, P[0], P[1], P[2], P[3]);
    flint_printf("pivot columns %wd %wd %wd, others %wd %wd\n",
                 piv[0], piv[1], piv[2], piv[3], piv[4]);
    nmod_mat_print_pretty(LU);

    /* unpack L (4 x r, unit diagonal) and U (r x 5) */
    nmod_mat_init(L, 4, r, 17);
    nmod_mat_init(U, r, 5, 17);
    for (i = 0; i < 4; i++)
        for (j = 0; j < r; j++)
            nmod_mat_entry(L, i, j) = (i == j) ? 1 :
                                      (j < i) ? nmod_mat_entry(LU, i, j) : 0;
    for (i = 0; i < r; i++)
        for (j = i; j < 5; j++)
            nmod_mat_entry(U, i, j) = nmod_mat_entry(LU, i, j);

    /* row i of L*U is row P[i] of A */
    nmod_mat_init(T, 4, 5, 17);
    nmod_mat_mul(T, L, U);
    nmod_mat_permute_rows(A, P, NULL);
    flint_printf("L*U == rows of A permuted by P: %d\n", nmod_mat_equal(T, A));

    nmod_mat_clear(A);
    nmod_mat_clear(LU);
    nmod_mat_clear(L);
    nmod_mat_clear(U);
    nmod_mat_clear(T);
    return 0;
}

rank 3, P = [1 0 2 3]
pivot columns 0 2 3, others 1 4
<4 x 5 integer matrix mod 17>
[ 2  4  1  0  5]
[ 0  0  3  1  4]
[ 9  5  0  2  4]
[ 1  1  0  0  0]
L*U == rows of A permuted by P: 1

Row 0 of $A$ has a zero in column 0, so rows 0 and 1 are swapped. Column 1 has no pivot because it is twice column 0. In row 2 of the output, the entries 9 and 5 belong to $L$ and the entries $0, 2, 4$ in columns 2 to 4 to $U$, whose pivot in this row is the 2 in column 3. Row 3 lies in the span of the rows above, so it holds only its row of $L$, which says that it is the sum of the first two rows of $U$.