Combining Algebraic and Numerical Techniques for Computing Matrix Determinant

Computing the sign of the determinant or the value of the determinant of an n × n matrix A is a classical well-know problem and it is a challenge for both numerical and algebraic methods. In this paper, we review, modify and combine various techniques of numerical linear algebra and rational algebraic computations (with no error) to achieve our main goal of decreasing the bit-precision for computing det A or its sign and enable us to obtain the solution with few arithmetic operations. In particular, we improved the precision     H p 2 log bits of the p -adic lifting algorithm ( H = 2 h for a natural number h ), which may exceed the computer precision β (see Section 5.2), to at most     β bits (see Section 6). The computational cost of the p -adic lifting can be performed in O ( hn 4 ). We reduced this cost to O ( n 3 ) by employing the faster p -adic lifting technique (see Section 5.3).


Introduction
Computation of the sign of the determinant of a matrix and even the determinant itself is a challenge for both numerical and algebraic methods.That is, to testing whether ( ) det 0 A > , ( ) det 0 A < , or det 0 A = for an n × n matrix A. In computational geometry, most decisions are based on the signs of the determinants.Among the geometric applications, in which the sign of the determinant needs to be evaluated, are the computations of con-vex hulls, Voronoi diagrams, testing whether the line intervals of a given family have nonempty common intersection, and finding the orientation of a high dimensional polyhedron.We refer the reader to [1]- [5], and to the bibliography therein for earlier work.These applications have motivated extensive algorithmic work on computing the value and particulary the sign of the determinant.One of the well-known numerical algorithms is Clarkson's algorithm.In [6], Clarkson uses an adaption of Gram-Schmidt procedure for computing an orthogonal basis and employs approximate arithmetics.His algorithm runs in time bits to represent the values for a d × d determinant with b-bit integer entries.On the other hand, the authors of [7] proposed a method for evaluating the signs of the determinant and bounding the precision in the case where 3 n ≤ .Their algorithm asymptotic worst case is worse than that of Clarkson.However, it is simpler and uses respectively b and ( ) 1 b + -bit arithmetic.In this paper we examine two different approaches of computing the determinant of a matrix or testing the sign of the determinant; the numerical and the algebraic approach.In particular, numerical algorithms for computing various factorizations of a matrix A which include the orthogonal factorization A QR = and the triangular factorizations A LU = , r A P LU = , and r c A P LUP = seem to be the most effective algorithms for computing the sign of ( ) det A provided that det A is large enough relative to the computer precision.Alternatively, the algebraic techniques for computing det A modulo an integer M based on the Chinese remainder theorem and the p-adic lifting for a prime p use lower precision or less arithmetic operations.In fact, the Chinese remainder theorem required less arithmetic operations than the p-adic lifting but with higher precision due to (9).We also demonstrate some effective approaches of combining algebraic and numerical techniques in order to decrease the precision of the computation of the p-adic lifting and to introduce an alternative technique to reduce the arithmetic operations.If det A is small, then obtaining the sign of the determinant with lower precision can be done by effective algebraic methods of Section 4. Although, the Chinese remainder algorithm approach requires low precision computations provided that the determinant is bounded from above by a fixed large number.This motivated us to generalize to the case of any input matrix A and to decrease the precision at the final stage of the Chinese remaindering at the price of a minimal increase in the arithmetic cost.Furthermore, in Section 5 we extend the work to an algorithm that computes ( ) det mod A M using low precision by relying on the p-adic lifting rather than the Chinese remaindering.

Definitions, Basic Facts, and Matrix Norms
Definition 1.For an n × n matrix A, the determinant of A is given by either ( ) ( ) ( ) matrix obtained by deleting the i-th row and j-th column of A. Fact 1. Well-known properties of the determinant include ( ) ( ) ( ) for any two matrices , If A is a triangular matrix of order n, then its determinant is the product of the entries on the main diagonal.
Definition 3.For a matrix A of size m × n, we define 1-norm: , and 2-norm:

Numerical Computations of Matrix Determinant
We find the following factorizations of an n × n matrix A whose entries are either integer, real, or rational by Gaussian elimination.
A LU = (1) . r c A P LUP = (3) We reduce A to its upper triangular form U. The matrix L is a unit lower triangular whose diagonal entries are all 1 and the remaining entries of L can be filled with the negatives of the multiples used in the elementary operations.P r is a permutation matrix results from row interchange and P c is a permutation matrix results from column interchange.
Fact 2. If P is a permutation matrix, then det 1 P = ± .Moreover, det 1 I = , where I is an identity matrix.

Triangular Factorizations
The following algorithm computes ( ) for Equation (3).One can easily output the sign of det A from the above algorithm.We only need to compute the sign of det r P , det L , detU , and det c P at stage 2 and multiply these signs.The value of det r P and det c P can be easily found by tracking the number of row or column interchanges in Gaussian elimination process which will be either 1 or −1 since those are permutation matrices.As L is a unit lower triangular matrix, det 1 L = .Finally, the computations of detU required 1 n − multiplications since U is an upper triangular matrix.Therefore, the overall arithmetic operations of Algorithm 1 will be dominated by the computational cost at stage 1, that is, ( )( ) multiplications, the same number for additions, subtractions, and comparisons, but

Orthogonal Factorization
, where Q T is the transpose of Q.
Lemma 1.If Q is an orthogonal matrix, then ( ) ( ) as a product of 1 n − matrices of Householder transformations, and the upper triangular matrix T where ( ) . Now for 0 1 t ≤ ≤ , we define the transformation ( ) ( ) = − .This proves that det , we have . The Givens rotations can also be used to compute the QR factorization, where is the total number of rotations, G j is the the j-th rotation matrix, and T Q A R = is an upper triangular matrix.Since the Givens algorithm computes Q as a product of the matrices of Givens rotations, then ( ) , , = for two fixed integers i and k, and for a pair of nonzero real numbers c and s that satisfies 2 2 1 c s If the Householder is used to find the QR, Algo- rithm 2 uses 1 n − multiplications at stage 3.It involves ( ) O k rows or columns of A are successively replaced by new ones.Therefore, in this case, the Givens algorithm is more effective., , , k m m m  be distinct pairwise relatively prime integers such that

Algebraic Computations of Matrix Determinant
, and let D denote an integer satisfying the inequality Then D is a unique integer satisfying ( 6) and (7).In addition, ( )( ) det mod A M based on the Chinese remainder theorem.)Input: An integer matrix A, an algorithm that computes ( ) satisfying (5) and are pairwise relatively prime, and 2) Compute the integers M i , s i , and y i as in (8) 1, , i k ∀ = .
The values of y i in (8) Suppose that the input matrix A is filled with integers and an upper bound det is known.Then we may choose an integer M such that and compute ( ) Hence, det A can be recovered by using (10).

Modular Determinant
( ) . In order to calculate ( ) , we choose a prime p that is bigger than 2 d and perform Gaussian elimination (2) on ( ) . This is the same as Gaussian elimination over  , except that when dividing by a pivot element a we have to calculate its inverse modulo p by the ex- tended Euclidean algorithm.This requires ( ) ( ) , where D is a diagonal matrix and R is a unit upper triangular matrix.The latest factorization can be computed by the Givens algorithm [9].If the entries of the matrix A are much smaller than p, then we do not need to reduce modulo p the results at the initial steps of Gaussian elimination.That is, the computation can be performed in exact rational arithmetics using lower precision.In this case, one may apply the algorithm of [10] and [11] to keep the precision low.The computations modulo a fixed integer 1 M > can be performed with precision  bits.In such a computation, we do not need to worry about the growth of the intermediate values anymore.However, to reduce the cost of the computations, one can work with small primes modular instead of a large prime, that is, these primes can be chosen very small with logarithmic length.Then the resulting algorithm will have lower cost of ( )

3
O n k and can be performed in parallel.Now, one can find a bound on the det A without actually calculating the determinant.This bound is given by Hadamard's inequality which says that ( ) where ∈  (nonnegative integers) is the maximal absolute value of an entry of A.

p-Adic Lifting of Matrix Inverses
In this section, we present an alternative approach for computing matrix determinant to one we discussed in earlier sections.The main goal is to decrease the precision of algebraic computations with no rounding errors.The technique relies on Newton's-Hensel's lifting (p-adic lifting) and uses ( ) n arithmetic operations.However, we will also show how to modify this technique to use order of n 3 arithmetic operations by employing the faster p-adic lifting.Given an initial approximation to the inverse, say S 0 , a well-known formula for Newton's iteration rapidly improves the initial approximation to the inverse of a nonsingular n × n matrix A:

S S I AS S S AS I S A S i
2) Finally, outputs h S .Note that, ( ) ∀ .This follows from the equation The above algorithm uses ( )

3
O n arithmetic operations performed modulo p J at the j-th stage, that is, a total of ( )

3
O hn arithmetic operations with precision of at most 2 log J p     bits.

p-Adic Lifting of Matrix Determinant
We can further extend p-adic lifting of matrix inverses to p-adic lifting of matrix determinants, that is ( ) det mod A p , using the following formula [12]: Here, A k denotes the k × k leading principal (north-western) submatrix of A so that n Moreover, ( ) , k k B denotes the ( ) , k k -th entry of a matrix B. In order to use Formula (14), we must have the inverses of A k available modulo a fixed prime integer p for all k .A nonsingular matrix A modulo p is called strongly nonsingular if the inverses of all submatrices A k exist modulo p.In general, a matrix A is called strongly nonsingular if A is nonsingular and all k × k leading principle submatrices are nonsingular.Here, we assume that A is strongly nonsingular for a choice of p. Finally, Algorithm 4 can be extended to lifting det A (see Algorithm 5).

Algorithm 5. (p-adic lifting of matrix determinant.)
Input: An integer 1 p > , an n × n matrix A, the matrices ( ) , and a natural number h.Output: ( )

Computations: 1) Apply Algorithm 4 to all pairs of matrices k
A and 0,k S , (replacing the matrices A and 0 S in the algo- rithm), so as to compute the matrices ( ) 3) Compute and output the value ( )( ) , as the reciprocal of ( ) The overall computational cost of Algorithm 5 at stage 1 is ( ) We apply Algorithm 4 (that is, ( ) , ) to all pairs of matrices: We then compute the value Therefore, we output the value of ( )

Faster p-Adic Lifting of the Determinant
Assuming A is strongly nonsingular modulo p, block Gauss-Jordan elimination can be applied to the block 2 × 2 matrix to obtain the well-known block factorization of A: The following is the divide-and-conquer algorithm for computing ( ) det A .Algorithm 6. (Computing the determinant of a strongly nonsingular matrix.)Input: An n × n strongly nonsingular matrix A. Output: 1) Partition A according to the block factorization above, where B is an  ) ( ) arithmetic operations using Algorithm 4. From the above factorization of A, we conclude that The computational cost of computing the determinant is ( ) ( ) . This can be derived from the above factorization of A using recursive factorization applied to B and S and the inverses modulo 2 H   p .Here, ( ) I k is the cost of computing the inverse and ( ) D k is the cost of computing the determinant.

Improving the Precision of the p-Adic Lifting Algorithm
Definition 6.Let e be a fixed integer such that where β is the computer precision.Then any integer z, such that 0 1 z e ≤ ≤ − , will fit the computer precision and will be called short integer or e-integer.The integers that exceed 1 e α − = will be called long integers and we will associate them with the e-polynomials ( ) , 0 i p e ≤ < for all i.Recall that Algorithm 4 uses 2 log j p     bits at stage j.For large j, this value is going to exceed β.In this section we decrease the precision of the p-adic algorithm and keep it below β.In order to do this, we choose the base where K is a power of 2, 2K H ≤ , K divides 2 H , and .Then for J K ≥ , we associate the p-adic lifting step of ( 13) with the matrix polynomial .
The input polynomial matrices are ( ) We perform the computations in (17) modulo J K x .The idea here is to apply e-reduction to all the entries of the output matrix polynomial followed by a new reduction modulo J K x .The resulting entries are just polynomials with integer coefficients ranging between 1 e − and 1 e − .This is due to the recursive e-reductions and then taking modular reductions again.Note that the output entries are polynomials with coefficients in the range γ − to γ for 2 β γ ≤ even before applying the e-reductions.This shows that the β-bit precision is sufficient in all of the computations due to (16).The same argument can be applied to the matrices , j k S for all j and k n < .

Numerical Implementation of Matrix Determinant
In this section we show numerical implementation of the determinant of n × n matrices based on the triangular factorization LU, P r LU, and P r LUP c .Algorithm 1 computes the triangular matrices , and let e′ , a, l  , and u  denote the maximum absolute value of the entries of the matrices E′ , A, L  , and U  .Also, E′ is the error matrix of the order of the roundoff in the entries of A.
Assuming floating point arithmetic with double precision (64-bits) and round to β -bit.Then the upper bound on the magnitude of the relative rounding errors is , where  is the machine epsilon.Our goal is to estimate e′ . .Then under (18),

(
) e e a nlu From the following matrix norm property , , , for 1, 2, q = ∞ , we obtain the bound .e E ∞ ′ ′ ≤ Theorem 3. Let a + denotes the maximum absolute value of all entries of the matrices ( ) i A computed by Gaussian elimination, which reduces A′ to the upper triangular form and let  as defined above.Then .

ln e e E n a an
By using the following results, we can estimate the magnitude of the perturbation of det A and det A′ caused by the perturbation of A′ .Here we assume that , We use the following facts to estimate d e of (21) in Lemma 5 below.
3) Substitute e + for e and ( )

2 .
where , i i r are the main diagonal entries of R. Here, we consider two well-known effective algorithms for computing the QR factorization; the House- holder and the Given algorithms.The Householder transformation (or reflection) is a matrix of the form H is symmetric and orthogonal.That is, + multiplications at stage 1, and ( ) O n evaluations of the square roots of positive numbers.If the Givens algorithm is used to find the QR, then it involves + multiplications at stage 1, and ( ) 2 O n evaluations.This shows some advantage of the Householder over the Givens algorithm.However, the Givens rotations algorithm allows us to update the QR of A successively by using only ( ) 2 O kn arithmetic operations and square root evaluations if ( )

By Algorithm 4 ,
compute the matrices:

refers to the floor of 2 n .) Compute 1 B − and S using the matrix equation 1 S
is strongly nonsingular modulo p, we can compute the inverse of k × k matrix by

− 1 J
with the e-polynomials in x for = 2 j J and for all j and k.These polynomials have degrees ( ) K − and take values of the entries for x e = .Define the polynomial matrices matrices r P  and c P  , where E L and E U are the perturbations of L and U.In general, we can assume that r

. 3 n
elimination with complete pivoting gives the matrix U rounded to 4 bits: and E L are the matrices obtained from accumulation of rounding errors.Also We now compute the upper bound e + of (20).Therefore, = is the size of the input matrix A .Hence we obtain the following upper bound d e + on d e by Lemma 5.That is, ( ) are computed by the Euclidean algorithm when applied to s i and m i for At stage 3, only one single multiplication is needed.All the above operations are calculated modulo p 2H .Now, we will estimate the value of H and p that must satisfy the bound22 det Apply Algorithm 1 in floating point arithmetic with unit roundoff bounded by  .Let triangular matrix, approximating the factor U of A from (2) or (3).