Computing determinants

fmpz_mat_det computes the determinant of a square integer matrix exactly. For small dimensions it uses cofactor expansion or fraction-free elimination, whose intermediate entries never exceed the size of minors of the matrix. For larger dimensions it computes the determinant modulo word-size primes, where elimination costs O(n^\omega) word operations per prime, and reconstructs it by the Chinese remainder theorem. The number of primes is governed by Hadamard's bound, unless a large divisor of the determinant is found first, in which case only the cofactor needs to be reconstructed. All results are proven unless proved is set to 0 in the modular functions.

Choice of algorithm#

For an n \times n matrix whose largest entry has b bits, fmpz_mat_det chooses as follows.

The empty matrix has determinant 1. A nonsquare matrix raises an exception in fmpz_mat_det; the other functions assume squareness.

Fraction-free elimination#

fmpz_mat_det_bareiss runs fmpz_mat_fflu with rank check on a copy of A. By Sylvester's identity the k-th pivot of Bareiss' algorithm is the leading k \times k minor of the row-permuted matrix, so the last pivot is \pm \det(A), with the sign of the permutation, and all intermediate entries are minors of A, bounded by Hadamard's inequality [Bar1968]. The cost is O(n^3) multiplications and exact divisions of integers of up to n(b + \log_2 n) bits, about O(n^5 b^2) word operations with classical integer multiplication, which is why it is used only for small n; see Fraction-free elimination.

Modular computation#

fmpz_mat_det_modular_given_divisor(det, A, d, proved) computes x = \det(A)/d for a known divisor d of the determinant (1 in fmpz_mat_det_modular):

1. H = fmpz_mat_det_bound(A) (Hadamard's bound), and the target is 2H/d, rounded up.
2. For primes p above 2^{59} that do not divide d: reduce A modulo p, compute its LU decomposition with nmod_mat_lu_with_pivots, and set x_p = \pm \prod_i U_{ii} \cdot d^{-1} \bmod p.
3. Combine the x_p with fmpz_CRT_ui in the symmetric range until the product of the primes exceeds 2H/d. Then |x| \le H/d is determined.

If the rank modulo some p is below n, the function tries to prove \det(A) = 0 at once with an exact kernel vector (_fmpz_mat_certify_singular_lu_mod_p), which succeeds whenever the rank modulo p equals the rank of A. With proved = 0 the loop also stops when the reconstructed value has stayed the same while more than 100 bits of primes were added: the probability that a wrong value survives that many primes is negligible for nonadversarial input, but the result is not proven.

The cost is about \log_2 (2H/d) / 59 LU decompositions modulo primes, each O(n^\omega), plus the reduction of A, O(n^2) reductions of b-bit entries per prime. With \log_2 H \approx n (b + \frac12 \log_2 n), this is O(n^{\omega + 1} (b + \log n)) word operations for fmpz_mat_det_modular.

Large divisors#

fmpz_mat_det_modular_accelerated first calls fmpz_mat_det_divisor, which solves Ax = v for a vector v with entries \pm 1 by p-adic lifting and returns the least common multiple d of the denominators of x. This d divides the largest invariant factor s_n of A, and for random matrices it equals s_n, which is \det(A) up to a small factor [AbbottBronsteinMulders1999; Pan1988]. Only \det(A)/d remains to be reconstructed, which typically takes one or two primes. Since the p-adic solution costs about the same as O(n) primes in the purely modular method but saves O(n(b + \log n))/59 primes, this pays off when the determinant is large compared with n, which fmpz_mat_det estimates by n \ge b. For matrices with a large Smith form cofactor \det(A)/s_n, for example diagonal-dominant structured matrices, the gain disappears; see Determinants.

Related functions#

• fmpz_mat_det_bound computes Hadamard's bound, and fmpz_mat_det_bound_submatrix a bound valid for all minors.
• fmpz_mat_is_singular decides \det(A) = 0 without computing the determinant, with a prime of 17 bits and an exact certificate; it is much faster than fmpz_mat_det for large nonsingular matrices.
• fmpz_mat_det_divisor returns the divisor described above.
• nmod_mat_det computes determinants modulo n, and fmpq_mat_det determinants of rational matrices by clearing denominators row by row.

Example#

A dense 40 \times 40 matrix with prescribed determinant 2^{40} \cdot 3^{20} \cdot 1000003, obtained from a sparse matrix by random unimodular row and column operations. All methods agree, and the Hadamard bound is far larger than the determinant. Since fmpz_mat_randdet spreads the prime factors of the determinant over many diagonal entries, the largest invariant factor, which is what fmpz_mat_det_divisor finds, is only 36 \cdot 1000003 here; for such matrices the accelerated method gains nothing.

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

int main(void)
{
    fmpz_mat_t A;
    fmpz_t D, d, t;
    flint_rand_t state;
    slong n = 40;

    flint_rand_init(state);
    fmpz_init(D);
    fmpz_init(d);
    fmpz_init(t);
    fmpz_mat_init(A, n, n);

    /* D = 2^40 * 3^20 * 1000003 */
    fmpz_ui_pow_ui(D, 3, 20);
    fmpz_mul_2exp(D, D, 40);
    fmpz_mul_ui(D, D, 1000003);

    fmpz_mat_randdet(A, state, D);
    fmpz_mat_randops(A, state, 2000);
    flint_printf("entries up to %wd bits\n", FLINT_ABS(fmpz_mat_max_bits(A)));

    fmpz_mat_det(d, A);
    flint_printf("det:                %{fmpz}\n", d);
    flint_printf("correct:            %d\n", fmpz_equal(d, D));

    fmpz_mat_det_bareiss(t, A);
    flint_printf("bareiss agrees:     %d\n", fmpz_equal(t, D));
    fmpz_mat_det_modular(t, A, 1);
    flint_printf("modular agrees:     %d\n", fmpz_equal(t, D));
    fmpz_mat_det_modular_accelerated(t, A, 1);
    flint_printf("accelerated agrees: %d\n", fmpz_equal(t, D));

    fmpz_mat_det_divisor(t, A);
    flint_printf("divisor:            %{fmpz}\n", t);
    flint_printf("divides det:        %d\n", fmpz_divisible(D, t));

    fmpz_mat_det_bound(t, A);
    flint_printf("bits of det %wd, of Hadamard bound %wd\n",
                 fmpz_bits(D), fmpz_bits(t));

    fmpz_mat_clear(A);
    fmpz_clear(D);
    fmpz_clear(d);
    fmpz_clear(t);
    flint_rand_clear(state);
    return 0;
}

entries up to 47 bits
det:                3833771493727452464601366528
correct:            1
bareiss agrees:     1
modular agrees:     1
accelerated agrees: 1
divisor:            36000108
divides det:        1
bits of det 92, of Hadamard bound 1865