Safe and Effective Determinant Evaluation

Kenneth L. Clarkson AT&T Bell Laboratories Murray Hill, New Jersey 07974

                e-mail:  clarkson@research.att.com

February 25, 1994

Abstract The problem of evaluating the sign of the determinant of a small matrix arises in many geometric algorithms. Given an n × n matrix A with integer entries, whose columns are all smaller than M in Euclidean norm, the algorithm given here evaluates the sign of the determinant det A exactly. The algorithm requires an arithmetic precision of less than 1.5n + 2 lg M bits. The number of arithmetic operations needed is O(n^{3}) + O(n^{2}) log OD(A)/β, where OD(A)| det A| is the product of the lengths of the columns of A, and β is the number of “extra” bits of precision,

min{lg(1/u) − 1.1n − 2 lg n − 2, lg N − lg M − 1.5n − 1},

where u is the roundoff error in approximate arithmetic, and N is the largest representable integer. Since OD(A) ≤ M^{n}, the algorithm requires O(n^{3} lg M) time, and O(n^{3}) time when β = Ω(log M).

1 Introduction

Many geometric algorithms require the evaluations of the determinants of small matrices, and testing the signs (±) of such determinants. Such testing is fundamental to algorithms for finding convex hulls, arrangements of lines and line segments and hyperplanes, Voronoi diagrams, and many others. By such numerical tests, a combinatorial structure is defined. If the tests are incorrect, the resulting structure may be wildly different from a “sensible” result[11, 6], and programs for computing the structures may crash.

Two basic approaches to this problem have been proposed: the use of exact arithmetic, and the design of algorithms to properly use inaccurate numerical tests. While the former solution is general, it can be quite slow: naively, n-tuple precision appears to be necessary for exact computation of the determinant of an n×n matrix with integer entries. While “adaptive” precision has been proposed and used[7], the resulting code and time complexity are still substantially larger than one might hope.

The second approach is the design (or redesign) of algorithms to use limited precision arithmetic and still guarantee sensible answers. (For example, [4, 2, 13, 6, 14, 12, 11, 9, 17].) Such an approach can be quite satisfactory, when obtainable, but seems applicable only on an ad hoc and limited basis: results for only a fairly restricted collection of algorithms are available so far.

This paper takes the approach of exact evaluation of the sign of the determinant, or more generally, the evaluation of the determinant with low relative error. However, the algorithm requires relatively low precision: less than 3n/2 bits, and additionally some number more bits than were used to specify the matrix entries. The algorithm is naturally “adaptive,” in the sense that the running time is proportional to the logarithm of the orthogonality defect OD(A) of the matrix A, where

OD(A) ≡ \frac{∏_{1≤i≤n}∥a_{i}∥}{| det A|}.

Here a_{i} is the ith column of A. (We let ∥a_{i}∥ ≡ \sqrt{a_{i}^{2}} denote the Euclidean norm of a_{i}, with a_{i}^{2} ≡ a_{i} · a_{i}.) In geometric problems, the orthogonality defect may be small for most matrices, so such adaptivity should be a significant advantage in running time. Note that the limited precision required by the algorithm implies that “native” machine arithmetic can be used and still allow the input precision to be reasonably high. For example, the algorithm can handle 10 × 10 matrices with 32-bit entries in the 53 bits available in double precision on many modern machines.

The algorithm is amenable to “rank-one updates,” so det B can be evaluated quickly (in about O(n^{2}) time) after det A is evaluated, if B is the same as A in all but one column.

One limitation of the approach here is that the matrix entries are required to be integers, rather than say rational numbers. This may entail preliminary scaling and rounding of input to a geometric algorithm; this limitation can be ameliorated by allowing input expressed in homogeneous coordinates, as noted by Fortune[1]. The resulting output is plainly the output for inputs that are “near” the originals, so such an algorithm is stable in Fortune’s sense.[2]

The new algorithm uses ideas from Lovász’s basis reduction scheme [8, 10, 16], and can be viewed as a “low-rent” version of that algorithm: a result of both algorithms is a new matrix A′, produced by elementary column operations, whose columns are roughly orthogonal. However, basis reduction does not allow a column of A to be replaced by an integer multiple of itself: the resulting set of matrix columns generates only a sublattice of the lattice generated by the columns of A. Since here only the determinant is needed, such a scaling of a column is acceptable since the scaling factor can be divided out of the determinant of the resulting matrix. Thus the problem solved here is easier than basis reduction, and the results are sharper with respect to running time and precision. The algorithm here yields a QR factorization of A, no matter what the condition of A, and so may be of interest in basis reduction, since Lovász’s algorithm starts by finding such a factorization. In practice, that initial factorization is found as an approximation using floating-point arithmetic; if the matrix MGS(matrix A) {

for k := 1 upto n do c_{k} := a_{k}; for j := k − 1 downto 1 do c_{k} −= a_{j}(c_{k}/c_{j});

return C; }

Figure 1: A version of the modified Gram-Schmidt procedure.

matrix is ill-conditioned, the factorization procedure may fail.

Schnorr has given an algorithm for basis reduction requiring O(n+lg M) bits; his algorithm for this harder problem is more complicated, and his constants are not as sharp.[15]

2 The algorithm

The algorithm is an adaptation of the modified Gram-Schmidt procedure for computing an orthogonal basis for the linear subspace spanned by a set of vectors. The input is the matrix A = [a_{1} a_{2} … a_{n}] with the column n-vectors a_{i}, i = 1 … n. One version of the modified Gram-Schmidt procedure is shown in Figure 1. We use the operation a/b for two vectors a and b, with a/b ≡ (a·b)/b^{2}. Note that (a + b)/c = a/c + b/c, and (a − (a/c)c) · c = 0. Implemented in exact arithmetic, this procedure yields an orthogonal matrix C = [c_{1} c_{2} … c_{n}], so c_{i} · c_{j} = 0 for i ≠ j, such that a_{k} = c_{k} + ∑_{1≤j<k} R_{jk}c_{j}, for some values R_{jk} = a_{k}/c_{j}. (In variance with some usage, in this paper orthogonal matrices have pairwise orthogonal columns: they need not have unit length. If they do, the matrix will be termed orthonormal.)

The vector c_{k} is initially a_{k}, and is then reduced by the c_{j} components of a_{k}: after step j of the inner loop, c_{k} ·c_{j} = 0. Note that even though a_{j} is used in the reduction, rather than c_{j}, the condition c_{k} · c_{i} = 0 for j < i < k is unchanged, since a_{j} has no c_{i} components for i > j.

Since A = CR where R is unit upper triangular, det A = det C; since C is an orthogonal matrix, | det C| = ∏_{1≤j≤n}∥c_{j}∥, and it remains only to determine the sign of the determinant of a perfectly conditioned matrix.

When the same algorithm is used with floating-point arithmetic to find a matrix B approximating C, the algorithm may fail if at some stage k a column b_{k} of B is very small: the vector a_{k} is very nearly a linear combination of the vectors {a_{1}, …, a_{k−1}}: it is close to the linear span of those vectors. Here the usual algorithm might simply halt and return the answer that the determinant is nearly zero. However, since b_{k} is small, we can multiply a_{k} by some large scaling factor s, reduce sa_{k} by its c_{j} components for j < k, and have a small resulting vector. Indeed, that vector will be small even if we reduce using rounded coefficients for a_{j}. With such integer coefficients, and using the columns float det safe(matrix A) {

float denom := 1; S^{b} := 0; for k := 1 upto n do

loop b_{k} := a_{k}; for j := k − 1 … 1 b_{k} −= a_{j}(b_{k}/b_{j});

if a_{k}^{2} ≤ 2b_{k}^{2} then {S^{b} += b_{k}^{2}; exit loop;}

s := S^{b}/4(b_{k}^{2} + 2δ_{k}a_{k}^{2});

s := min\{s, (N/1.58^{k} − \sqrt{S^{b}})^{2}/a_{k}^{2}\}; if s < 1/2 then s := 1 else s := ⌈\sqrt{s} ⌋; if s = 1 and a_{k}^{2} ≤ 0.9S^{b} then s := 2;

denom ∗= s; a_{k} ∗= s; for j := k − 1 … 1 a_{k} −= a_{j}⌈(a_{k}/b_{j})⌋; if a_{k} = 0 then return 0;

end loop; return det approx(B)/denom;

}

Figure 2: An algorithm for estimating det A with small relative error.

a_{j} and exact arithmetic for the reduction steps, det A remains unchanged, except for being multiplied by the scaling factor s.

Thus the algorithm given here proceeds in n stages as in modified Gram-Schmidt, processing a_{k} at stage k. First the vector b_{k} is found, estimating the component of a_{k} orthogonal to the linear span of {a_{1}, …, a_{k−1}}. A scaling factor s inversely proportional to ∥b_{k}∥ is computed, and a_{k} is replaced by sa_{k} and then reduced using exact elementary column operations that approximate Gram-Schmidt reduction. The result is a new column a_{k} = sĉ_{k} + ∑_{1≤j<k} α_{j}c_{j}, where ĉ_{k} is the old value of that vector, and the coefficients α_{j} ≈ 1/2. Hence the new a_{k} is “more orthogonal” to earlier vectors than the old one. The processing of a_{k} repeats until a_{k} is nearly orthogonal to the previous vectors a_{1}, …, a_{k−1}, as indicated by the condition a_{k}^{2} ≤ 2b_{k}^{2}.

The algorithm is shown in Figure 2. All arithmetic is approximate, with unit roundoff u, except the all-integer operations a_{k} ∗= s and the reduction steps for a_{k} (although a_{k}/b_{j} is computed approximately). The value N is no more than the largest representable integer. For a real number x, the integer ⌈x⌉ is the least integer no smaller than x, and ⌈x⌋ is ⌈x − 1/2⌉. The quantity δ_{k} is defined in Notation 3.4 below.

The precision requirements of the algorithm can be reduced somewhat by “reorthogonalization,” that is, simply applying the reduction step to b_{k} again after computing it, or just before exiting the loop for stage k. Hoffman discusses and analyzes this technique; his analysis can be sharpened in our situation.[5]

Note that the condition a_{k}^{2} ≤ 2b_{k}^{2} may hold frequently or all the time; if the latter, the algorithm is very little more than Modified Gram-Schmidt, plus the procedure for finding the determinant of the very-well-conditioned matrix B.

3 Analysis

The analysis requires some elementary facts about the error in computing the dot product and Euclidean norm.

Lemma 3.1 For n-vectors a and b and using arithmetic with roundoff u, a · b can be estimated with error at most 1.01nu∥a∥∥b∥, for nu ≤ .01. The relative error in computing a^{2} is therefore 1.01nu, under these conditions. The termination condition for stage j implies a_{j}^{2} ≤ 2b_{j}^{2}(1 + 2.03nu).

Proof: For a proof of the first statement, see [3], p. 35. The remaining statements are simple corollaries.

We’ll need a lemma that will help bound the error due to reductions using b_{j}, rather than c_{j}.

Lemma 3.2 For vectors a, b, and c, let d ≡ b − c, and let δ ≡ ∥d∥/∥b∥. Then

a/b − a/c ≤ ∥a∥∥d∥/∥c∥∥b∥ = δ∥a∥/∥c∥.

Note also a/b ≤ ∥a∥/∥b∥.

Proof: The Cauchy-Schwartz inequality implies the second statement, and also implies

a/b − a/c ≤ ∥a∥∥b/b^{2} − c/c^{2}∥.

Using elementary manipulations

\frac{b}{b^{2}} − \frac{c}{c^{2}} = \frac{(dc^{2} − cd · c) − c(d · c + d^{2})}{c^{2}b^{2}}.

The two terms in the numerator are orthogonal, so the norm of this vector is

\frac{\sqrt{(dc^{2} − cd · c)^{2} + c^{2}(d · c + d^{2})^{2}}}{c^{2}b^{2}} = \frac{\sqrt{c^{2}d^{2}b^{2}}}{c^{2}b^{2}} = δ/∥c∥,

as desired.

Lemma 3.3 For vectors a, b and c, and d, and δ as in the previous lemma, suppose a^{2} ≤ 2b^{2}(1 + α), a/c = 1, and α < 1. Then ∥a − b∥/∥b∥ ≤ 1 + 4δ + α.

Proof: We have

\frac{(a − b)^{2}}{b^{2}} = \frac{a^{2} − 2a · b + b^{2}}{b^{2}} ≤ 3 + 2α − 2a · b/b^{2},

using a^{2} ≤ 2b^{2}(1 + α). By the assumption a/c = 1,

\frac{a · b}{b^{2}} = \frac{a · (c + d)}{b^{2}} = \frac{c^{2} + a · d}{b^{2}} ≥ (1 − δ)^{2} − \sqrt{2(1 + α)},

where the last inequality follows using the triangle and Cauchy-Schwartz inequalities and the assumption a^{2} ≤ 2b^{2}(1 + α). Hence

\frac{(a − b)^{2}}{b^{2}} ≤ 3 + 2α − 2((1 − δ)^{2} − \sqrt{2(1 + α)}δ) ≤ 1 + 2α + 8δ,

giving ∥a − b∥/∥b∥ ≤ 1 + 4δ + α, as claimed.

Notation 3.4 Let

η ≡ 3(n + 2)u,

and for j = 1 … n, let

δ_{j} ≡ L_{j}η − η,

where

L_{j}^{2} ≡ jμ^{2}(μ^{2} + 2.04)^{j−1},

and

μ ≡ 1.58.

Let

ϕ ≡ 1 + 2.04/μ^{2}.

Let

d_{k} ≡ b_{k} − c_{k}

for k = 1 … n.

Theorem 3.5 Suppose the reductions for b_{k} are performed using rounded arithmetic with unit roundoff u. Assume that (n + 2)u < .01 and δ_{n} < 1/32. (For the latter, lg(1/u) ≥ 1.1n + 2 lg n + 7 suffices.) Then:

(i) ∥b_{k} − a_{j}(b_{k}/b_{j})∥ ≤ 1.54∥b_{k}∥.

(ii) The computed value of a single reduction step for b_{k}, or b_{k} − (b_{k}/b_{j})a_{j}, differs from the exact value by a vector of norm no more than η∥b_{k}∥. Before a reduction of b_{k} by a_{j}, ∥b_{k}∥ ≤ ∥a_{k}∥μ^{k−1−j}.

(ii) After the reduction loop for b_{k}, d_{k} satisfies

∥d_{k}∥ ≤ δ_{k}∥a_{k}∥/\sqrt{2},

and after stage k, ∥d_{k}∥ ≤ δ_{k}∥b_{k}∥.

Proof: First, part (i): let

a′_{j} ≡ (a_{j} − b_{j}) − ((a_{j} − b_{j})/b_{j})b_{j},

which is orthogonal to b_{j}, and express b_{k} − a_{j}(b_{k}/b_{j}) as

b_{k} − (b_{k}/b_{j})b_{j}
− (b_{k}/b_{j})a′_{j}
− (b_{k}/b_{j})((a_{j} − b_{j})/b_{j})b_{j}.

Since b_{k}^{2} = α^{2} + β^{2}b_{j}^{2} + γ^{2}(a′_{j})^{2}, for some α, β = b_{k}/b_{j}, and γ = b_{k}/a′_{j}, we have

(b_{k} − (b_{k}/b_{j})(b_{j} + a′_{j}))^{2} = α^{2} + (γ − β)^{2}(a′_{j})^{2},

and so

(b_{k} − (b_{k}/b_{j})(b_{j} + a′_{j}))^{2}/b_{k}^{2} ≤ (γ − β)^{2}/(β^{2}b_{j}^{2}/(a′_{j})^{2} + γ^{2}),

which has a maximum value

1 + (a′_{j})^{2}/b_{j}^{2} ≤ 2 + 4δ_{j} + 2.03nu,

and so

∥b_{k} − (b_{k}/b_{j})(b_{j} + a′_{j})∥ ≤ 1.47∥b_{k}∥.  (1)

Using part (ii) inductively and Lemmas 3.3 and 3.1, and the assumption δ_{j} ≤ 1/32,

∥(b_{k}/b_{j})((a_{j} − b_{j})/b_{j})b_{j}∥

≤ \frac{∥b_{k}∥}{b_{j}^{2}}(a_{j} − c_{j} − d_{j}) · (c_{j} + d_{j})

= \frac{∥b_{k}∥}{b_{j}^{2}}(−c_{j} · d_{j} + d_{j} · (a_{j} − b_{j}))

≤ ∥b_{k}∥(δ_{j}(1 + δ_{j}) + δ_{j}(1 + 4δ_{j} + 2.03nu))

≤ ∥b_{k}∥δ_{j}(2.15)

Thus with (1),

∥b_{k} − a_{j}(b_{k}/b_{j})∥ ≤ (1.47 + 2.15/32)∥b_{k}∥ ≤ 1.54∥b_{k}∥.

This completes part (i).

Now for part (ii). From Lemma 3.1 and Lemma 3.2, the error in computing b_{k}/b_{j} = b_{k} · b_{j}/b_{j}^{2} is at most

((1 + u)(1 + 2.03nu) − 1)∥b_{k}∥/∥b_{j}∥,

and ∥a_{j}∥ ≤ \sqrt{2}(1+1.02nu)∥b_{j}∥, and so the difference between a_{j}(b_{k} ·b_{j}/b_{j}^{2}) and its computed value is a vector ϵ with norm no more than ∥b_{k}∥ times

\sqrt{2}(1 + 1.02nu)((1 + u)^{2}(1 + 2.03nu) − 1)  (2)
< 2.9002(n + 1)u  (3)

A reduction step computes the nearest floating-point vector to b_{k} −a_{j}(b_{k}/b_{j})+ϵ. Hence the error is a vector ϵ + ϵ′, where

∥ϵ′∥ ≤ u∥b_{k} − a_{j}(b_{k}/b_{j}) + ϵ∥
≤ u(∥b_{k} − a_{j}(b_{k}/b_{j})∥ + ∥ϵ∥).

With this and (2), the computed value of b_{k} after the reduction step, differs from its real value by a vector that has norm no more than ∥b_{k}∥ times

2.9002(n + 1)u(1 + u) + 1.54u
< 3(n + 2)u
= η.

Using part (i), the norm of b_{k} after the reduction step is no more than (1.54 + η)∥b_{k}∥ < μ∥b_{k}∥.

We turn to part (iii), using (ii) and inductively (iii) for j < k. (We have d_{1} = 0 for the inductive basis.)

The vector b_{k} is initially a_{k}, which is c_{k} plus a vector in the linear span of a_{1} … a_{k−1}, which is the linear span of c_{1} … c_{k−1}. When b_{k} is reduced by a_{j}, it is replaced by a vector

b_{k} − (b_{k}/b_{j})a_{j} + ϵ_{j},  (4)

where ϵ_{j} is the roundoff. Except for ϵ_{j}, the vector b_{k} − c_{k} continues to be in the linear span of {c_{1} … c_{k−1}}. Hence d_{k} ≡ b_{k} − c_{k} can be expressed as e_{k} + Υ, where e_{k} = ∑_{1≤j<k} γ_{j}c_{j}, for some values γ_{j}, and Υ is no longer than the sum of the roundoff vectors, so

∥Υ∥ ≤ η∥a_{k}∥μ^{k−1}/(μ − 1),  (5)

using part (ii). From the triangle inequality,

d_{k}^{2} ≤ (∥e_{k}∥ + Υ)^{2}.  (6)

The quantity γ_{j} is

\frac{1}{c_{j}} · [b_{k} − (b_{k}/c_{j})a_{j} + (b_{k}/c_{j})a_{j} − (b_{k}/b_{j})a_{j}]

= b_{k}/c_{j} − b_{k}/b_{j},

where b_{k} is the value of that vector for the a_{j} reduction. (Note that except for roundoff, b_{k}/c_{j} does not change after the reduction by a_{j}.) By Lemma 3.2 and this part of Theorem 3.5 as applied to d_{j},

|γ_{j}| ≤ ∥b_{k}∥δ_{j}/∥c_{j}∥,  (7)

Using part (ii) and (7),

e_{k}^{2} = ∑_{1≤j<k} γ_{j}^{2}c_{j}^{2}
≤ ∑_{1≤j<k} δ_{j}^{2}μ^{2(k−1−j)}a_{k}^{2}.

Using (5),(6), and δ_{j} < L_{j}η,

\frac{d_{k}^{2}}{a_{k}^{2}} ≤ μ^{2k} \frac{η^{2}}{μ^{2}} (\frac{1}{μ − 1} + \sqrt{∑_{1≤j<k} L_{j}^{2}/μ^{2j}})^{2} .  (8)

It suffices for part (iii) that

δ_{k} ≥ μ^{k}\sqrt{2(1 + η)}\frac{η}{μ} (\frac{1}{μ − 1} + \sqrt{∑_{1≤j<k} L_{j}^{2}/μ^{2j}}) ,

where the (1 + η) term allows the statement

d_{k}^{2} ≤ δ_{k}^{2}a_{k}^{2}/2(1 + 2.03nu).  (9)

Letting M_{k} = L_{k}/μ^{k}, we need

M_{k} ≥ \frac{1}{μ^{k}} + \frac{\sqrt{2(1 + η)}}{μ} (\frac{1}{μ − 1} + \sqrt{∑_{1≤j<k} M_{j}^{2}}) ,

Since L_{1} = 1 allows

∥d_{1}∥ = 0 = δ_{1} = η − η,

M_{k} = \sqrt{j(1 + 2(1 + η)/μ^{2})^{k−1}}, or L_{k}^{2} = kμ^{2}(μ^{2} + 2.04)^{k−1} for k > 1, gives sufficiently large δ_{k}, using the facts

∑_{1≤j<k} jx^{j−1} = \frac{kx^{k−1}}{x − 1} − \frac{x^{k} − 1}{(x − 1)^{2}}

and

\sqrt{x − y} ≤ \sqrt{x} − y/2\sqrt{x}.

The last statement of part (iii) follows from the termination condition for stage k, which from Lemma 3.1 implies a_{k}^{2} ≤ 2b_{k}^{2}(1 + 2.03nu) < 2b_{k}^{2}(1 + η).

Lemma 3.6 With Notation 3.4 and the assumptions of Theorem 3.5, let S_{k}^{c} ≡ ∑_{1≤j<k} c_{j}^{2} and S_{k}^{b} the computed value of ∑_{1≤j<k} b_{j}^{2}. Then

|S_{k}^{c} − S_{k}^{b}|/S_{c}^{k} < 1/10.

Proof: With Theorem 3.5(iii), we have

b_{j}^{2}/c_{j}^{2} ≤ 1/(1 − δ_{j})^{2},

implying

∑_{1≤j<k} b_{j}^{2}/S_{k}^{c} ≤ 1/(1 − δ_{j})^{2}.

Lemma 3.1 implies S_{k}^{b}/ ∑_{1≤j<k} b_{j}^{2} ≤ 1 + 1.03(k + n)u, and so

\frac{S_{k}^{b}}{S_{k}^{c}} ≤ \frac{1 + 1.03(k + n)u}{(1 − δ_{j})^{2}},

implying S_{k}^{b} − S_{k}^{c} < S_{k}^{c}/10. A similar argument shows that S_{k}^{c} − S_{k}^{b} < S_{k}^{c}/10. For the remaining discussion we’ll need even more notation:

Notation 3.7 Let

f_{k} ≡ a_{k}%c_{k} ≡ a_{k} − (a_{k}/c_{k})c_{k},

so that a_{k} = c_{k} + f_{k}, c_{k} and f_{k} are perpendicular, and a_{k}^{2} = c_{k}^{2} + f_{k}^{2}. Suppose a reduction loop for a_{k} is done. Let â_{k} denote a_{k} before scaling it by s, and let ĉ_{k} and f̂_{k} denote the components c_{k} and f_{k} at that time. Let a_{k}^{(j)} denote the vector a_{k} when reducing by a_{j}, so a_{k}^{(k−1)} is a_{k} before the reduction loop, and a_{k}^{(0)} is a_{k} after the reduction loop. We’ll use z_{jk} to refer to a computed value of a_{k}/b_{j}.

Lemma 3.8 With assumptions as in the previous theorem, before a reduction of a_{k} by a_{j},

∥a_{k}∥ ≤ s∥â_{k}∥μ^{k−1−j} + (3/4) ∑_{j<i<k}∥c_{i}∥μ^{i−1−j}.

Proof: When reducing by a_{j}, a_{k} = a_{k}^{(j)} is replaced by

a_{k}^{(j−1)} = a_{k} − z_{jk}a_{j} + αa_{j},

for some |α| ≤ 1/2. As in Theorem 3.5(ii), the norm of a_{k} − z_{jk}a_{j} is no more than μ∥a_{k}∥, so

∥a_{k}^{(j−1)}∥ ≤ μ∥a_{k}^{(j)}∥ + ∥a_{j}∥/2.

Thus

∥a_{k}^{(j)}∥ ≤ s∥â_{k}∥μ^{k−1−j} + ∑_{j<i<k}∥a_{i}∥μ^{i−1−j}/2
≤ s∥â_{k}∥μ^{k−1−j} + (3/4) ∑_{j<i<k}∥c_{i}∥μ^{i−1−j},  (10)

where the last inequality follows from Lemma 3.1 and Theorem 3.5(ii).

Lemma 3.9 With the assumptions of Theorem 3.5, in a reduction loop for a_{k},

ĉ_{k}^{2} ≤ 0.55â_{k}^{2}

and

f̂_{k}^{2} ≥ 0.45â_{k}^{2}.

Proof: Since the reduction loop is done, the conditional “a_{k}^{2} ≤ 2b_{k}^{2}” returned false. Therefore, using Theorem 3.5(iii) and reasoning as in Lemma 3.1,

∥â_{k}∥ > \sqrt{2}∥b_{k}∥(1 − 2.03nu)
≥ \sqrt{2}(∥ĉ_{k}∥ − δ_{k}∥â_{k}∥/\sqrt{2})(1 − 2.03nu),

which implies

ĉ_{k}^{2} ≤ 0.55â_{k}^{2},  (11)

using the assumption δ_{k} < 1/32. The last statement of the lemma follows from â_{k}^{2} = ĉ_{k}^{2} + f̂_{k}^{2}.

Theorem 3.10 With assumptions as Theorem 3.5, after the reduction loop for a_{k} it satisfies

a_{k}^{2} ≤ max{â_{k}^{2}, 2.6S_{k}^{c}}.

Also

S_{k}^{c} ≤ M^{2}(3.3)^{k−2}.

The condition

lg N ≥ lg M + 1.5n + 1

allows the algorithm to execute.

Proof: We begin by bounding the c_{j} components of a_{k} for j < k, and use these to bound ∥f_{k}∥, after the reduction loop. Note that the c_{j} component of a_{k} remains fixed after the reduction of a_{k} by a_{j}. That is, using Notation 3.7,

a_{k}^{(0)}/c_{j} = a_{k}^{(j−1)}/c_{j} = (a_{k}^{(j)} − ⌈z_{jk}⌋a_{j})/c_{j}.

Using a_{j}/c_{j} = 1, this implies

|a_{k}^{(0)}/c_{j}|
≤ |a_{k}^{(j)}/c_{j} − a_{k}^{(j)}/b_{j}| + |a_{k}^{(j)}/b_{j} − ⌈z_{jk}⌋|
≤ ∥a_{k}^{(j)}∥δ_{j}/∥c_{j}∥ + (1/2 + 2.2(n + 1)u∥a_{k}^{(j)}∥/∥c_{j}∥)
< 1/2 + ∥a_{k}^{(j)}∥(δ_{j} + η)/∥c_{j}∥,  (12)

where the next-to-last inequality follows from Theorem 3.5(iii) and from reasoning as for (2). Since

f_{k}^{2} = ∥a_{k}^{(0)}∥^{2} − s^{2}ĉ_{k}^{2} = ∑_{1≤j<k}(a_{k}^{(0)}/c_{j})^{2}c_{j}^{2},

and the triangle inequality implies ∥x+y∥^{2} ≤ (∥x∥+∥y∥)^{2}, we have f_{k}^{2} no more than

∑_{1≤j<k}(∥c_{j}∥/2 + ∥a_{k}^{(j)}∥(δ_{j} + η))^{2}

≤ (\sqrt{S_{k}^{c}}/2 + \sqrt{∑_{1≤j<k}(∥a_{k}^{(j)}∥(δ_{j} + η))^{2}})^{2}.  (13)

From the definition of δ_{j} and (10),

∑_{1≤j<k}∥a_{k}^{(j)}∥^{2}(δ_{j} + η)^{2}

≤ η^{2}k ∑_{1≤j<k} (s∥â_{k}∥μ^{k−j} + \frac{3}{4} ∑_{j<i<k}∥c_{i}∥μ^{i−j})^{2}(μ^{2}ϕ)^{j}

≤ η^{2}k (\sqrt{∑_{1≤j<k} s^{2}â_{k}^{2}μ^{2k}ϕ^{j}}

+ \frac{3}{4}\sqrt{∑_{1≤j<k} ( ∑_{j<i<k}∥c_{i}∥μ^{i})^{2}ϕ^{j}})^{2}
≤ δ_{k}^{2}(s∥â_{k}∥μ/\sqrt{2} + 3\sqrt{S_{k}^{c}}/4)^{2}

With (13), f_{k}^{2} is no more than

(\sqrt{S_{k}^{c}}/2 + δ_{k}(s∥â_{k}∥μ/\sqrt{2} + 3\sqrt{S_{k}^{c}}/4))^{2}
≤ (\sqrt{S_{k}^{c}}(1/2 + δ_{k}) + 1.2δ_{k}s∥â_{k}∥)^{2} .  (14)

If s = 1, then using Lemma 3.9 and δ_{k} ≤ 1/32,

a_{k}^{2} = c_{k}^{2} + f_{k}^{2} ≤ 0.55â_{k}^{2} + (\sqrt{S_{k}^{c}}(1/2 + δ_{k}) + 1.2δ_{k}∥â_{k}∥)^{2}
≤ 0.55â_{k}^{2} + (0.54\sqrt{S_{k}^{c}} + 1.2∥â_{k}/32∥)^{2}
≤ max{â_{k}^{2}, S_{k}^{c}}.  (15)

If s = 2 because the conditional “s = 1 and a_{k}^{2} ≤ 0.9S^{b}” returned true, then ∥â_{k}∥ ≤ 1.01\sqrt{S_{k}^{c}}, and

a_{k}^{2} ≤ 4(0.55)â_{k}^{2} + (\sqrt{S_{k}^{c}}(1/2 + δ_{k}) + 2(1.2)δ_{k}∥â_{k}∥)^{2}  (16)
≤ 2.6S_{k}^{c}.  (17)

Now suppose s > 1 and the conditional “s = 1 and a_{k}^{2} ≤ 0.9S^{b}” returned false. By Theorem 3.5(iii), ĉ_{k}^{2} ≤ b_{k}^{2} + 2δ_{k}â_{k}^{2}, so that when ⌈\sqrt{s}⌋ is evaluated it is no more than \sqrt{S^{b}/4c_{k}^{2}}(1.03) + 1/2, where the 1.03 factor bounds the error in computing s. Thus

c_{k}^{2} ≤ s^{2}ĉ_{k}^{2} ≤ (\sqrt{S^{b}/4}(1.03) + ∥ĉ_{k}∥/2)^{2}
≤ (0.55\sqrt{S_{k}^{c}} + ∥ĉ_{k}∥/2)^{2}
≤ (0.55\sqrt{S_{k}^{c}} + 0.375∥â_{k}∥)^{2},  (18)

using Lemma 3.9. When ⌈\sqrt{s}⌋ is evaluated, it is smaller than

1.03\sqrt{S^{b}/4δ_{k}â_{k}^{2}} + 1/2,

and so

sδ_{k}∥â_{k}∥ ≤ 0.55\sqrt{δ_{k}S_{k}^{c}} + δ_{k}∥â_{k}∥/2.

With (14) and the assumption δ_{k} ≤ 1/32,

f_{k}^{2} ≤ (\sqrt{S_{k}^{c}}(1/2 + δ_{k}) + 0.55\sqrt{δ_{k}S_{k}^{c}} + δ_{k}∥â_{k}∥/2)^{2}
≤ (0.63\sqrt{S_{k}^{c}} + δ_{k}∥â_{k}∥/2)^{2},  (19)

and so with (18),

a_{k}^{2} = c_{k}^{2} + f_{k}^{2}
≤ (0.55\sqrt{S_{k}^{c}} + 0.375∥â_{k}∥)^{2} + (0.63\sqrt{S_{k}^{c}} + ∥â_{k}∥/64)^{2}
≤ max{â_{k}^{2}, 1.4S_{k}^{c}}.  (20)

With (15), (16), and (20), we have a_{k}^{2} ≤ max{â_{k}^{2}, 2.6S_{k}^{c}}. First we show that c_{k}^{2} ≤ max{ĉ_{k}^{2}, 2.2S_{k}^{c}}. It remains to bound S_{k}^{c} inductively. If s = 1 then c_{k}^{2} = ĉ_{k}^{2}. If s = 2 because the conditional “s = 1 and a_{k}^{2} ≤ 0.9S^{b}” returned true, then c_{k}^{2} ≤ 2.2S_{k}^{c}. If s > 1 and the conditional “s = 1 and a_{k}^{2} ≤ 0.9S^{b}” returned false, then with (18),

c_{k}^{2} ≤ (0.55\sqrt{S_{k}^{c}} + ∥ĉ_{k}∥/2)^{2} ≤ max\{ĉ_{k}^{2}, 2.2S_{k}^{c}\},

Let σ_{k}^{c} be an upper bound for S_{k}^{c}, so σ_{2}^{c} = M^{2} suffices, and when stage k is finished, σ_{k+1}^{c} = σ_{k}^{c} + 2.2σ_{k}^{c} is large enough, so σ_{k}^{c} = (3.2)^{k−2}M^{2} is an upper bound for S_{k}^{c}.

Theorem 3.11 Let β be the minimum of

lg(1/u) − 1.1n − 2 lg n − 2

and

lg N − lg M − 1.5n − 1.

(The conditions of Theorems 3.5 and 3.10 are implied by β > 4.) The algorithm requires O(n^{3}) + O(n^{2}) log OD(A)/β time.

Proof: Clearly each iteration of the main loop requires O(n^{2}) time. The idea is to show that for each stage k, the loop body is executed O(1) times except for executions that reduce OD(A) by a large factor; such reductions may occur when a_{k} is multiplied by a large scale factor s, increasing det A, or by reducing ∥a_{k}∥ substantially when s is small.

We consider some cases, remembering that c_{k}^{2} never decreases during stage k. The discussion below assumes that S^{b}/4(b_{k}^{2} + 2δ_{k}a_{k}^{2}) evaluates to no more than (N/1.58^{k} − \sqrt{S^{b}})^{2}/a_{k}^{2}; if not, s is limited by the latter, and β is bounded to the second term in its definition.

s > 2, b_{k}^{2} ≥ δ_{k}â_{k}^{2}. After two iterations with this condition, either s ≤ 2 or c_{k}^{2} will be sufficiently large that 2b_{k}^{2} ≥ a_{k}^{2}.

s > 2, b_{k}^{2} ≤ δ_{k}â_{k}^{2}. Here (19) and Lemma 3.9 show that

∥f_{k}∥ ≤ 0.62\sqrt{S_{k}^{c}} + δ_{k}∥â_{k}∥/2,

and s > 2 implies

δ_{k}∥â_{k}∥/2 ≤ δ_{k}\sqrt{S^{b}/δ_{k}}/\sqrt{4(2.5)2}/2 < 0.04\sqrt{S_{k}^{c}},

so ∥â_{k}∥ ≤ 0.66\sqrt{S_{k}^{c}}/\sqrt{.45} < \sqrt{S_{k}^{c}} during the second iteration under these conditions. Thus s ≥ 0.37/\sqrt{δ_{k}} here for all but possibly the first iteration, while a_{k}^{2} is bounded. Hence lg OD(A) decreases by 0.5 lg(1/δ_{k}) − 2 during these steps.

s ≤ 2. From (14) and Lemma 3.9,

∥f_{k}∥ ≤ (1/2 + δ_{k})\sqrt{S_{k}^{c}} + 2.4δ_{k}∥f̂_{k}∥/\sqrt{.45}
≤ 0.54\sqrt{S_{k}^{c}} + 4delta_{k}∥f̂_{k}∥,

so y_{k} ≡ ∥f_{k}∥/0.54\sqrt{S_{k}^{c}} converges to 1/(1 − 4δ_{k}).

Suppose y_{k} ≥ 1/(4δ_{k})^{2}/(1 − 4δ_{k}); here s = 1, and y_{k} decreases by a factor of 4δ_{k}/(1 − 4δ_{k}) ≤ 4.6δ_{k} at each reduction loop. Hence ∥a_{k}∥ is reduced by a factor of 4.6δ_{k}, while det A remains the same. Hence lg OD(A) decreases by at least lg(1/δ_{k}) − 2.3.

Suppose y_{k} < 1/(4δ_{k})^{2}/(1 − 4δ_{k}); then in O(1) (less than 7) reductions, y_{k} ≤ \sqrt{1.01}/(1 − 4δ_{k}), so that

a_{k}^{2} ≤ f_{k}^{2}/0.45
≤ 1.01(0.54)^{2}S_{k}^{c}/(1 − 4δ_{k})^{2}/0.45
≤ 0.75S_{k}^{c},

and so the conditional “a_{k}^{2} ≤ 0.9S^{b}” will return true. Thereafter O(1) executions suffice before stage k completes, as c_{k} doubles in magnitude at each step.

Here is one way to get the estimate det approx(B): estimate the norms ∥b_{i}∥, normalize each column of B to obtain a matrix B̄ with columns b̄_{j} ≡ b_{j}/∥b_{j}∥, for j = 1, …, n, and then return det B̄ ≈ ±1 times ∏_{1≤i≤n}∥b_{i}∥. To suggest that det B̄ can be estimated using low-precision arithmetic, we show that indeed | det B̄| ≈ 1 and that the condition number of B̄ is small. In the proof, we use the matrix C̄ whose columns are c̄_{j} ≡ c_{j}/∥c_{j}∥, for j = 1, …, n. (Note that the truly orthonormal matrix C̄ has det C̄ = 1 and Euclidean condition number κ_{2}(C̄) = 1.)

Theorem 3.12 With assumptions as in Theorem 3.5, for any x ∈ R^{n},

|∥B̄^{T} x∥/∥x∥ − 1| ≤ 1.68δ_{n},

and so | det B̄| ≥ 1 − 2nδ_{n}, and B̄ has condition number κ_{2}(B̄) ≤ 1 + 3δ_{n}.

Proof: It’s easy to show that for x ∈ R^{n},

|b̄_{j} · x − c̄_{j} · x| ≤ (δ_{j} + η)∥x∥,

where the η on the right accounts for roundoff in computing b̄_{j} from b_{j}. Thus

(B̄^{T} x − C̄^{T} x)^{2} ≤ x^{2} ∑_{1≤j≤n}(δ_{j} + η)^{2}
≤ x^{2}η^{2}n ∑_{1≤j≤n}(μ^{2} + 2.04)^{j−1}
≤ 2δ_{n}^{2}x^{2}.

Since ∥C̄^{T} x∥ = ∥x∥,

|∥B̄^{T} x∥/∥x∥ − 1| ≤ \sqrt{2}δ_{n}.

Recall that B̄^{T} and B̄ have the same singular values.[3] Since the bound above holds for any x, the singular values σ_{j} of B̄ satisfy |σ_{j} −1| ≤ δ_{n}, for j = 1, …, n, and so

κ_{2}(B̄) = σ_{1}/σ_{n} ≤ \frac{1 + \sqrt{2}δ_{n}}{1 − \sqrt{2}δ_{n}} ≤ 1 + 3δ_{n},

using (1 + x)/(1 − x) = 1 + 2x/(1 − x) and 1/(1 − x) increasing in x. Since | det B̄| = σ_{1}σ_{2} ⋯ σ_{n}, we have

| det B̄| ≥ (1 − \sqrt{2}δ_{n})^{n} ≥ exp(−1.5nδ_{n}) ≥ 1 − 1.5nδ_{n},

where we use 1 − x ≥ exp(−x\frac{log(1−α)}{−α}) for 0 < x ≤ α < 1.

With these results, it is easy to show that Gaussian elimination with partial pivoting can be used to find det B̄, using for example the backwards error bound that the computed factors L and U have LU = B̄ +Δ, where Δ is a matrix with norm no more than n^{2}\sqrt{n}2^{n−1}∥B̄∥u; this implies, as above, that the singular values, and hence the determinant, of LU are close to B.

Acknowledgements

It’s a pleasure to thank Steve Fortune, John Hobby, Andrew Odlyzko, and Margaret Wright for many helpful discussions.

References

[1] S. Fortune. personal communication.

[2] S. Fortune. Stable maintenance of point-set triangulations in two dimensions. In Proc. 30th IEEE Symp. on Foundations of Computer Science, pages 494–499, 1989.

[3] G. H. Golub and C. F. Van Loan. Matrix Computations. Johns Hopkins University Press, Baltimore and London, 1989.

[4] D. Greene and F. Yao. Finite-resolution computational geometry. In Proc. 27th IEEE Symp. on Foundations of Computer Science, pages 143–152, 1986.

[5] W. Hoffman. Iterative algorithms for Gram-Schmidt orthogonalization. Computing, 41:335–348, 1989.

[6] C. Hoffmann, J. Hopcroft, and M. Karasick. Robust set operations on polyhedral solids. IEEE Comp. Graph. Appl., 9:50–59, 1989.

[7] M. Karasick, D. Lieber, and L. Nackman. Efficient delaunay triangulation using rational arithmetic. ACM Transactions on Graphics, 10:71–91, 1990.

[8] A. K. Lenstra, H. W. Lenstra, and L. Lovász. Factoring polynomials with rational coefficients. Math. Ann., 261:515–534, 1982.

[9] Z. Li and V. Milenkovic. Constructing strongly convex hulls using exact or rounded arithmetic. In Proc. Sixth ACM Symp. on Comp. Geometry, pages 235–243, 1990.

[10] L. Lovász. An algorithmic theory of numbers, graphs, and complexity. SIAM, Philadelphia, 1986.

[11] V. Milenkovic. Verifiable implementations of geometric algorithms using finite precision arithmetic. Artificial Intelligence, 37:377–401, 1988.

[12] V. Milenkovic. Verifiable Implementations of Geometric Algorithms using Finite Precision Arithmetic. PhD thesis, Carnegie Mellon U., 1988.

[13] V. Milenkovic. Double precision geometry: a general technique for calculating line and segment intersections using rounded arithmetic. In Proc. 30th IEEE Symp. on Foundations of Computer Science, pages 500–505, 1989.

[14] D. Salesin, J. Stolfi, and L. Guibas. Epsilon geometry: building robust algorithms from imprecise calculations. In Proc. Seventh ACM Symp. on Comp. Geometry, pages 208–217, 1989.

[15] C. P. Schnorr. A more efficient algorithm for lattice basis reduction. J. Algorithms, 9:47–62, 1988.

[16] A. Schrijver. Theory of Linear and Integer Programming. Wiley, New York, 1986.

[17] K. Sugihara and M. Iri. Geometric algorithms in finite-precision arithmetic. Technical Report RMI 88-10, U. of Tokyo, 1988.