Perhaps the most important class of expressions in geometric computation is the determinant. We examine some algorithmic techniques needed for computing determinants in the context of exact geometric computation.

The main issue is to compute some approximation to the determinant which has the correct sign. Clarkson [4] gave a general method for n×n matrices with L-bit integer entries, assuming limited precision exact integer arithmetic. The method can be very efficient when combined with machine floating point arithmetic for n and L that are moderate. When n = 2 or 3, Avnaim et al [1] gave efficient solutions when the inputs are machine representable. In [2], Brönnimann and Yvinec revisit Clarkson's algorithm and generalizes the method of Avnaim et al.

Filters for Determinants.   If we use a standard Gaussian elimination algorithm for computing the determinant, then a simple filter is just to compare the absolute value of the computed determinant with the well-known error bounds found in Numerical Analysis books. Pan and Yanquiang [7] pointed out that this filter (which is static) is not very efficacious. They propose instead to use the the distance from the nearest singular instance of a matrix. This requires some extra but not too expensive computation, but leads to much better filtering ratios.

1   Concepts from Numerical Analysis

In this section, we recall some basic concepts developed in numerical analysis. We assume vectors in \Bbb Cⁿ and matrices in \Bbb C^m×n for various natural numbers n,m ³ 1. For a matrix A Î \Bbb C^m×n, then A^* Î \Bbb C^n×m denotes its Hermitian transpose. In the special case where A is real, A^* is equal to A^T, the transpose. Let |A| be the matrix whose (i,j)th entry is |a_ij|. Two real matrices A, B of the same dimension can be compared componentwise: we write A £ B if a_ij £ b_ij for all i,j. The identity matrix is denoted I (of appropriate dimension determined by context).

References: John [6], Higham [5], Chaitin-Chatelin and Frayssé [3].

Vector Norms.   We want some measure of ``size'' of a vector x Î \Bbb Cⁿ. This is needed, for instance, in estimating computational errors. This is captured by the notion of a norm. A norm on \Bbb Cⁿ is a function N:\Bbb Cⁿ® \Bbb R such that for all x,y Î \Bbb Cⁿ,

• N(x) ³ 0, with equality iff x = 0.
• N(ax) = |a| N(x) for all a Î \Bbb C.
• N(x+y) £ N(x)+N(y)

||x||p : = Ö[p] |x1|p + ¼+|xn|p

where x = (x₁ ,¼, x_n)^T. The special cases of p = 1 and p = 2 are noteworthy. In particular, the 2-norm is differentiable and invariant under an orthogonal transformation of the space. We can also generalize this to p = ¥ and define

||x||¥ : = max { |x1| ,¼, |xn| }.

Two norms N and N¢ are equivalent if there exists positive constants c < C such that for all x Î \Bbb Cⁿ,

c·N(x) £ N¢(x) £ C·N(x).

It is easy to see that this notion of equivalence is transitive.

Lemma 1 Any two vector norms are equivalent.

Proof.  By transitivity of equivalence, it suffices to show that any norm N is equivalent to the ¥-norm. Write x = å_{i = 1}ⁿ x_i eⁱ where eⁱ is the ith elementary vector. Then

N(x) £ n å i = 1  |xi| N(ei) £ ||x||¥ n å i = 1  N(ei) £ C ||x||¥

where C = å_{i = 1}ⁿ N(eⁱ). Now, consider the unit sphere under the ¥-norm, S = {x Î \Bbb Cⁿ: ||x||_¥ = 1}. Since S is compact, the norm function N:S® \Bbb R achieves its mininum value at some x⁰ Î S. Let N(x⁰) = c. If ||x||_¥ = b then we have

N(x) ³ N(b · x b ) = b N( x b ) ³ b·c = c ||x||¥.

We have thus specified the requisite constants c,C in the definition of equivalence. Q.E.D.

Matrix Norms.   We can treat a matrix A Î \Bbb C^m×n as a mn-vector and apply the previous vector norms. In particular, the vector 2-norm when applied to A in this way is called the Frobenious norm, and given the special notation ||A||_F. If we view A as a linear operator from \Bbb Cⁿ to \Bbb C^m, then we can use a general method for defining norms on linear operators. To do this, we fix two vector norms N^m and Nⁿ on \Bbb C^m and \Bbb Cⁿ, respectively. We define the derived norm N¢ on A Î \Bbb C^m×n as follows:

N¢(A) : = sup x ¹ 0, x Î \Bbb Cn  Nm(Ax) Nn(x) .

(1)

Since N(lx) = |l|·N(x), we can equivalently define N¢ as

N¢(A) : = sup Nn(x) = 1  Nm(Ax).

(2)

Thus N¢ measures the maximum ``expansion'' of a vector on the unit sphere {x Î \Bbb Cⁿ: Nⁿ(x) = 1} in \Bbb Cⁿ. We say that N¢ is derived from N^m, Nⁿ. The following properties are easy to show:

1. N¢(A) ³ 0 with equality iff A = \0.
2. N¢(aA) = |a|·N¢(A) for a Î \Bbb C.
3. N¢(A+B) £ N¢(A)+N¢(B).
4. N¢(AB) £ N¢(A)N¢(B).
5. N^m(Ax) £ N¢(A) Nⁿ(x).

Assume that Nⁿ and N^m are vector p-norms. We want to characterize the derived matrix norm, denoted N_p, for p = 1,2,¥.

Theorem 2 Let A = (a_ij) Î \Bbb C^m×n. Assume B : =
A^*A has eigenvalues m₁, m₂ ,¼, m_n.

(a) N¥(A) = max i  å j  |aij| (b) N1(A) = max j  å i  |aij| (c) N2(A) =   æ Ö max i  mi

Proof.  Let x = (x₁ ,¼, x_n)^T Î \Bbb Cⁿ be a non-zero vector.

(a) We have

||Ax||¥ = max i  | å j  aijxj| £ ||x||¥ max i  å j  |aij|.

The inequality becomes an equality if we choose x_j = |a_i0j|/a_i0j where i₀ is the index that achieves the maximum in max_i å_j |a_ij|. This proves N_¥(A) = max_i å_j |a_ij|. (b) Again, we have

||Ax||1 = å i  | å j  aijxj| £ (1) å i  å j  |aijxj| £ (2) å j  |xj| å i  |aijxj| £ (3) ||x||1 max j  å i  |aij|.

Suppose j₀ is the index that maximizes the expression max_j å_i |a_ij|, and we choose x = e^j₀ (the unit vector with 1 in the j₀th position). Then inequalities £ ₍₁₎ and £ ₍₂₎ becomes equalities because x is now a unit vector, while inequality £ ₍₂₎ becomes an equality because of our choice of j₀. This proves N₁(A) = max_j å_i |a_ij|. (c) For any m and x, if Bx = mx then mx^* x = x^*Bx = (x^*A^*)(Ax). Hence

m = ||Ax||22 ||x||22 .

This shows that the eigenvalues m_i of B are real and non-negative. It is well-known that for any Hermitian matrix B (i.e., B^* = B), there exists unitary a matrix T (i.e., T^*T = I) such that B = T^* D T where D is the diagonal matrix Diag(m₁ ,¼, m_n). Without loss of generality, assume m₁ is the largest eigenvalue. For any x where ||x||₂ = 1, we have

||Ax||22 = (Ax)*(Ax) = x* B x = x*T* D T x = y* D y       (y = Tx) = å i  mi |yi|2 £ m1       (||x||2 = ||y||2 = 1).

The final inequality is an equality if y = e¹ (the unit vector with 1 in the first component). This proves N₂(A) = Ö{max_i m_i}. Q.E.D.

We read (a) and (b) as saying: N_¥(A) is the maximum of the 1-norm of the rows of A; N₁(A) is the maximum of the 1-norm of the columns of A. To remember whether it is row or column, use this device: the symbol ``¥'' has a dominant horizontal direction, suggesting rows; likewise, the symbol ``1'' suggests columns.

It is not easy to compute N₂(A). Hence, the following connection between N₂(A) and the Frobenious norm ||A||_F of A is useful.

Lemma 3

Ön N2(A) ³ ||A||F ³ N2(A).

Proof.  Note that

||A||F2 = å i  å j  |aij|2 = Tr(A*A).

Writing B = A^*A = T^*DT as in the above proof, we see

Tr(B) = å i  mi.

Therefore

n max i  mi ³ Tr(B) ³ max i  mi.

The result follows by taking square-roots. Q.E.D.

These bounds on N₂(A) can be computed in O(n²) arithmetic operations.

2   Elements of Error Analysis

Let u be the unit roundoff for the floating point number system under discussion (see the chapter on floating point numbers). In the floating point system F(b, t), u = b^1-t/2. Recall that the standard error model says that the approximate value [x\tilde] of any floating point arithmetic operation is related to the exact value x by the relation

~ x   = x(1+d),        |d| £ u.

We develop some basic tools for the error analysis within this model. The following is a slight extension of Higham [5,p. 69].

In the following, we assume that n is small enough so that

nu < 1.

(3)

For any n ³ 1 and real e, define

gn(e) : = ne 1-ne .

(4)

When e = u, we simply write

gn : = gn(u).

Note that

1+gn(e) = 1 1-ne

(5)

Lemma 4  
(i) g_n(e) is a strictly increasing function for e < 1/n.

(ii) If 1 > 2nu then |g_n(e)| £ g_n iff e Î [ [(-u)/(1-2nu)],u].

(iii) Assume 1 > 2nu. For any real q, we have |q| £ g_n iff q = g_n(e) for some e Î [ [(-u)/(1-2nu)],u].

Proof.  (i) If e < e¢ < 1/n, then we verify that 1+g_n(e) = [1/(1-ne)] < [1/(1-ne¢)] = 1+g_n(e¢). Note that e, e¢ can be negative.

(ii) Since 1 > 2nu, we have [(-u)/(1-2nu)] < 0 < u. We check that g_n([(-u)/(1-2nu)]) = -g_n(u) = -g_n. By part (i), it follows that for all e Î [ [(-u)/(1-2nu)],u], we have -g_n £ g_n(e) £ g_n.

(iii) This is a consequence of (ii). Q.E.D.

Lemma 5 Let |d_i| £ u, r_i = ±1 and 1 > 2nu. Then

n Õ i = 1  (1+di)ri = 1 + gn(e)

where for some e such that |g_n(e)| £ u.

Proof.  Suppose n = 1. Then

(1+d1)r1 = ì ï ï í ï ï î 1+d1 if  r1 = 1, 1+ -d1 1+d1 if  r1 = -1 = ì ï ï í ï ï î 1+ g1 æ ç è d1 1+d1 ö ÷ ø if  r1 = 1, 1+ g1(-d1) if  r1 = -1 = 1+ g1(e)

where e Î [ [(-u)/(1-2u)],u]. Thus |g₁(e)| £ u, as required. Let n ³ 2. By induction,

n-1 å i = 1  (1+di)ri = 1+gn-1(e)

and

(1+dn)rn = 1+g1(e¢)

where |g_n-1(e)| £ u and |g₁(e¢)| £ u. It suffices to show that

1+gn(e¢¢) = (1+gn-1(e))(1+g1(e¢))

(6)

for some e¢¢ satisfying |g_n(e¢¢)| £ u. Alternatively, we must show that e¢¢ lies in the range [-u/(1-2nu), u], given that e Î [-u/(1-2(n-1)u),u] and e¢ Î [-u/(1-2u),u]. It follows from (5) that (6) amounts to showing

1 1-ne¢¢ = 1 1-(n-1)e · 1 1-e¢ 1-ne¢¢ = (1-(n-1)e) (1-e¢) e¢¢ = n-1 n e(1-e¢) + 1 n e¢.

Thus e¢¢ is a convex combination of e(1-e¢) and e¢. There are three cases to consider:

Case A. If e¢ ³ 0, then both e(1-e¢) and e¢ lie in the range [-u/(1-2nu), u]. So its convex combination e¢¢ lies in the range [-u/(1-2nu), u].

Case B. If e¢ < 0 and e ³ 0 then

e¢¢ £ n-1 n e(1-e¢) £ n-1 n u æ ç è 1- u 1-2u ö ÷ ø £ u

and

e¢¢ ³ 1 n e¢ ³ -u n(1-2u) ³ -u 1-2nu .

Case C. Finally, assume e¢ < 0 and e < 0. Then e¢¢ < 0 and we need to show that e¢¢ ³ -u/(1-2nu). Since e¢¢ is a convex combination of e(1-e¢) and e¢, it suffices to show that both e(1-e¢) and e¢ are at least -u/(1-2nu). Clearly, e¢ ³ -u/(1-2u) ³ -u/(1-2nu). Also,

e(1-e¢) ³ e(1+u) ³ -u(1+u) 1-2(n-1)u ³ -u 1-2nu ,

where the last inequality may be verified directly. Q.E.D.

On Evaluation of Expressions.   We consider expressions E in the basic arithmetic operations, generalized to allow summation (å_{i = 1}ⁿ) or multi-product (Õ_{i = 1}ⁿ). E.g.,

E = 3 x2 + x3 - (x1 x2 x3) + 5 å i = 1  xi.

An expression gives a mathematical value, but we also are interested in its evaluation on a computer. In this case, the order of evaluating the expression is important. This order is unique for a fully parenthesized expression. But in general, for any E, we define the natural evaluation of E to be the direct left-to-right of expressions, respecting any parentheses when present. Each summation or multiproduct notation is linearized in the obvious way, by increasing indices and evaluated left-to-right. E.g., å_{i = 1}ⁿ x_n is simply a shorthand for x₁+x₂+¼+x_n, whose natural evaluation order is

((¼((x1+x2)+ x3) + ¼+ xn-1)+xn)

Suppose E¢ is any equivalent expression, i.e., obtained from E by applying any valid laws of arithmetic (associativity, commutativity, distributivity, etc). Because of finite precision arithmetic, the evaluation of E¢ and E need not produce the same result. The simplest example is

E = x+y+z,        E¢ = x+z+y.

(7)

Another example is the following, using the law of distributivity:

E = (x+y)(v+w),        E¢ = xv + xw + yv + yw.

(8)

In numerical analysis, we tend not to exploit distributivity because they can render E¢ very different from E. Comparing the errors in the natural evaluations E and E¢ can be complicated. One hint of this is the fact that the number of operations in E and E¢ can be exponentially different:

E = Õ (1+xi),       E¢ = 1 + å i  xi + å i < j  xixj +¼ Õ i  xi.

Therefore, we focus on one class of mathematically equivalent expressions, called re-ordering equivalent. Two expressions are re-ordering equivalent if they can be obtained from each other by applying the laws of associativity and commutativity only. If E is re-ordering equivalent to E¢, then their natural evaluation incurs the same number of arithetic operations.

An approximate value [E\tilde] for E is called a natural value (resp. re-ordering value) if is obtained by the natural evaluation (resp.  re-ordering evaluation) of E.

Lemma 6 Let E = å_{i = 1}ⁿ x_i y_i.
(i) The natural value [E\tilde] of E satisfies [E\tilde] = å_{i = 1}ⁿ x_i y_i(1+q_i) where |q₁| £ g_n and for i = 2 ,¼, n, |q_i| £ g_n-i+2.
(ii) If E¢ is any re-ordering equivalent expression, its natural value [(E¢)\tilde] satisfies [(E¢)\tilde] = å_{i = 1}ⁿ x_i y_i(1+q_i) where each |q_i| £ g_n.

Proof.  Let E_j = å_{i = 1}^j x_i y_i for j = 1 ,¼, n.
(i) If n = 1, [E\tilde] = x₁ y₁(1+q₁), |q₁| £ g₁. This has the required form. For n ³ 2, we inductively have

~ En-1   = n-1 å i = 1  xiyi(1+q¢i)

where |q¢₁| £ g_n-1 and for i = 2 ,¼, n-1, |q¢_i| £ g_n-i+1. From

~ En   = ( ~ En-1   + xn yn (1+q¢n))(1+q) = n å i = 1  xiyi(1+q¢i)(1+q) = n å i = 1  xiyi(1+qi)

where |q_i| satisfies the lemma.
(ii) This is immediate from part (i). Q.E.D.

Lemma 7 Let E = (z - å_{i = 1}^n-1 x_i y_i)/y_n. If E¢ is any re-ordering equivalent expression, its natural value [(E¢)\tilde] satisfies

~ E¢   (1+qn) = z - n-1 å i = 1  xi yi(1+qi) yn

where each |q_i| £ g_n.

Proof.  Q.E.D.

3   Conditioning

Let f: X® Y be a continuous function on the normed spaces X, Y. The condition number k_f(x) of an input instance x Î X for problem f measures the sensitivity of f(x) to perturbations in x. In particular, ``ill or well conditioning'' refers to the relative largeness or smallness of k_f(x). We can work with an absolute or a relative notion of conditioning. We follow [8] in the following definitions.

There is an absolute and a relative concept of conditioning. The relative notion is more important, but since the absolute notion is simpler, we also will present it. The absolute conditioning function of f is defined by

kAf(x) : = lim d® 0  sup ||dx|| £ d  ||df|| ||dx|| .

(9)

Here, we write df for f(x+dx)-f(x). If f is differentiable, write J(x) = J_f(x) for the Jacobian of f at x. The (i,j)th entry of J(x) is ¶f_i/ ¶x_j(x). Then the quantities dx and df becomes (respectively) the infinitesimals dx and J(x)dx. Thus

kAf(x) = || Jf(x)||

(10)

where || J_f(x)|| is the derived norm for the norms on X and Y.

The relative conditioning function of f is defined as

kfR(x) : = lim d® 0  sup ||dx|| £ d  æ ç è ||df|| || f(x) || / ||dx|| ||x|| . ö ÷ ø

(11)

Again when f is differentiable, we get

kfR(x) = || Jf(x)|| || x|| || f(x)|| .

(12)

In the following, we simply write k_f(x) to refer to the relative conditioning function.

Condition Number of a Matrix.   Let f: \Bbb Cⁿ ® \Bbb C^m where

f(x) = Ax

where A is a m×n matrix. Applying the definitions directly, we obtain

kA(x) : = lim d® 0  sup ||dx|| £ d  æ ç è ||df|| || f(x)|| / ||dx|| ||x|| ö ÷ ø = lim d® 0  sup ||dx|| £ d  æ ç è ||A(x + dx) - Ax|| ||Ax|| / ||dx|| ||x|| ö ÷ ø = lim d® 0  sup ||dx|| £ d  æ ç è ||A dx|| ·||x|| ||Ax|| ·||dx|| ö ÷ ø = ||A|| ·||x|| ||Ax||

If A is square, and A^-1 exists, then we have || x || = ||A^-1Ax || £ ||A^-1|| ·||Ax||, Plugging for ||x|| into the above condition number, we obtain

kA(x) £ || A || ·|| A-1 || .

The nice thing about the righthand side is that it is independent of x. This quantity is so useful that it deserves its own name and notation:

k(A) : = || A || ·|| A-1 ||

is the condition number of the square matrix A.

This number k(A) is also when we want to solve the linear equation Ax = b where b is given and we need to compute x. This is just the original problem with A^-1 in place of A, but because of the way we define condition numbers of a matrix, the condition number of A can serve to bound the conditioning of this problem as well.

Condition Number of matrix. This is defined to be

k(A) : = ||A-1|| ||A|| .

Normwise vs. componentwise norm. |A| = [ |a_ij| ]

Distance to a Singularity.   The numerical stability of a numerical problem is influenced by its distance to the nearest singularity. For linear systems, this distance can be defined as follows: for a nonsingular matrix A and any matrix norm N, the closest distance from A to a singular matrix is defined as

d(A) : = inf S  {N(S-A)}

where S range over all singular matrices. We next characterize this distance. This result, attributed to Turing, was known to Banach in the 1920s. Its generalization to arbitrary norms is from Gastinel (1966). First, we note a simple fact:

Lemma 8 For any square matrix A and any derived norm N, if N(A) < 1 then I+A is invertible.

Proof.  This is equivalent to showing that (I+A)x = \0 has the unique solution x = \0. If (I+A)x = \0 then x = -Ax and hence N(x) = N(Ax) £ N(A)N(x). We claim that N(x) = 0 since otherwise N(A) < 1 implies N(x) £ N(A)N(x) < N(x). Thus x = \0. Q.E.D.

Theorem 9 [Turing-Gastinel] For non-singular A and derived norm N,

d(A) = inf S  { N(A-S)} = 1 N(A-1) .

Proof.  We first prove that d(A) ³ [1/(N(A^-1))]. Let S be any matrix such that N(S-A) < [1/(N(A^-1))]. It suffices to that S is non-singular. Writing D for S-A, we have N(A^-1D) £ N(A^-1)N(D) < 1. Thus (I+A^-1D) is invertible and we see that

A(I+A-1D) = A+D = S

As a product to two non-singular matrices, S is non-singular.

To show that d_T(A) £ [1/(N(A^-1)N(A) )], we construct a singular matrix S such that N(S-A) = [1/(N(A^-1))]. We pick two vectors x, y such that

y* A-1 x = N(A-1),       N(xy*) = 1.

We check that

æ ç è A - xy* N(A-1) ö ÷ ø A-1x = 0.

This proves that S = A - [(xy^*)/(N(A^-1))] is singular. Moreover, N(A-S) = N(xy^*)/N(A^-1) = 1/N(A^-1). Q.E.D.

Component-wise Distance.   For a non-singular square matrix A, define (see [3])

dT(A) : = inf S  ì í î N(S-A) N(A) ü ý þ

(13)

where S ranges all singular matrices. Thus d_T(A) is the relative distance from A to the nearest singular matrix S. The subscript refers to Turing. Rewrite the definition of d_A(T) in (13) as

dT(A) = inf S  { a: N(S-A) = aN(A) }

where S ranges over all singular matrices. We can further rewrite it as

dT(A) = inf S,a  { a: N(S-A) £ aN(A) }

The componentwise analogue of this definition is therefore

dC(A) : = inf S, a  { a: |S-A| £ a|A| }.

A characterization of d_C(A) is given by Rohn (1990) (see [3,p.66]).

4   Convex Hulls and Partial Compilation of Expressions

Sometimes, part of the arguments of an expression are fixed while others are vary. We can exploit this by constructing ``partially compiled expressions''. To see an example of this phenomenon, we look at the Beneath-Beyond Algorithm for convex hulls.

See [Page 147ff in Edelsbrunner].

Let S be a set of points in d dimensions. FOR SIMPLICITY, ASSUME GENERAL POSITION. We sort the points by their first coordinate, and add them in this order. The method (called Beneath-Beyond) begins with a d+1 simplex. At each step, we add a new point p to the current hull H:

Prove the 5 properties stated for the norm N¢ derived from N^m, Nⁿ.

Let A be a square matrix and N a matrix norm.
(i) Show that N(A) ³ |m| for any eigenvalue of A.
(ii) Let us show that the above inequality cannot be improved to an equality. Let

A = 1 1 0 1 .

Then the eigenvalues of A satisfy |m₁| = |m₂| = 1. Prove that for any norm N, we have N(A) > 1. Hint: Write A = I+E where I is the identity matrix. What is A^k for any k Î \Bbb N? Show that if N(A) £ 1 then N(E) = 0, contradiction.
(iii) (Fritz John) Show that for all square matrices A and all e > 0, there exists a norm N such that N(A) £ e+ max_m |m| where m ranges over the eigenvalues of A.

(a) Give specific numerical values for the expressions E and E¢ (7) so that their natural evaluations result is different answers. Assume IEEE arithmetic.

(b) Do the same for (8).

Improve lemma 5 in case r_i = 1 for all i.

(Trefethan) Determine the absolute and relative conditioning function of the following functions.

(i) f(x) = x/2 for x Î \Bbb C.  
(ii) Suppose f(x) = Öx for x > 0.  
(ii) Suppose f(x) = x-y for (x,y) Î \Bbb C². Use the ¥-norm on \Bbb C².

5  OLD

REFERENCE: Certified Computation of the Sign of a Matrix Determinant, Dario A. Bini, Victor Y. Pan and Yanqiang Yu. 1998, abstract.

An approach to computing the sign of a determinant D of an n×n matrix M is to numerically compute an approximate value A for D, estimate the roundoff error E, and to certify the computed sign of A in case |A| > E. Unfortunately, |D| varies from 0 to ||M||ⁿ, and the the known techniques for error analysis gives error bounds that are large even when |D| is small. A predecessor to this paper by Pan, Yu and Stewart uses this approach. The present paper avoids estimating the roundoff error E, but instead solve the simpler task of estimating the mininum distance 1/||M^-1|| from M to a singular matrix. The method uses O(n³) single precision flops. Based on this estimate, in case the sign is not certified, a decision can be made for increasing the precision of the numerical computation or to fall back to symbolic guaranteed methods. Thus, this algorithm can be used as a kind of filter.

Other work: Clarkson proposed a general algorithm based on Gram-Schmidt orthogonalization. The certification method is rather complicated to program. For n £ 4, Avnaim et al proposed an efficient ``single-precision'' method. Modular methods have been proposed by Brönnimann et al.

6  Roundoff Errors in Gaussian Elimination

Let A be a n×n matrix, and A_k be its kth column and a_i,j its (i,j)th entry. The following matrix operator norms are used: the 2-norm ||A||₂, the row norm ||A||_¥ = max_i å_j |a_i,j|, the column norm ||A||₁ = max_j å_i |a_i,j|. Also, let |A| denote the matrix with (i,j)th entry |a_i,j|. A triangular matrix is unit (resp., proper) if its diagonal entries are all 1 (resp., 0).

Theorem 10 (Higham, Theorem 9.3, page. 175): Let

A = PA¢Q,    ~ A   = A¢+ E = ~ L   ~ U

where the matrices are all n square, P and Q are permutation matrices, [L\tilde] and [U\tilde] are unit lower triangular and U upper triangular. If [L\tilde] and [U\tilde] are computed numerically from A by Gaussian elimination with complete pivoting, then

|E| £ gn | ~ L   |·| ~ U   |

where g_n = ne/(1-ne) and e is the unit roundoff.

Since det(A) = det(P)det(A¢)det(Q) and det(P) = ±1 and det(Q) = ±1 are known, it is enough to determine the sign of det(A¢).

Corollary 11 For any operator norm,

||E|| £ e = || (|| ~ L   ||·|| ~ U   ||) || gn

Computing the columns and row norms uses O(n²) arithmetic operations.

References

[1]

F. Avnaim, J.-D. Boissonnat, O. Devillers, F. Preparata, and M. Yvinec. Evaluating signs of determinants using single-precision arithmetic. Algorithmica, 17(2):111-132, 1997.

[2]

H. Brönnimann and M. Yvinec. Efficient exact evaluation of signs of determinants. Algorithmica, 27:21-56, 2000.

[3]

F. Chaitin-Chatelin and V. Frayssé. Lectures on Finite Precision Computations. Society for Industrial and Applied Mathematics, Philadelphia, 1996.

[4]

K. L. Clarkson. Safe and effective determinant evaluation. IEEE Foundations of Computer Science, 33:387-395, 1992.

[5]

N. J. Higham. Accuracy and stability of numerical algorithms. Society for Industrial and Applied Mathematics, Philadelphia, 1996.

[6]

F. John. Lectures on Advanced Numerical Analysis. Gordon and Breach, 150, Fifth Avenue, New York, NY 10011, 1967.

[7]

V. Y. Pan and Y. Yu. Certified computation of the sign of a matrix determinant. Proc. 10th ACM-SIAM Symp. on Discrete Algorithms (SODA99), pages 715-724, 1999. Accepted for Algorithmica.

[8]

L. N. Trefethen and I. David Bau. Numerical Linear Algebra. Society for Industrial and Applied Mathematics, Philadelphia, 1997.

Footnotes:

¹ We rely on the fact that these properties fit familiar ``mathematical patterns'' such the triangular inequality.