§1. Sylvester’s Identity  Lecture X  Page 1

Lecture X Linear Systems

Determinants and solving linear systems of equations are fundamental in algebraic computation. The basic facts of this topic are well-known when the underlying algebraic structure R is a field. Gaussian elimination is the main method in this case. But our main interest is when R is a domain D. Of course, we can still use embed a determinant computation in the quotient field Q_{D} of D. But this method turns out to be inefficient, for reasons similar to those that argue against computing the GCD over D[X] via the Euclidean algorithm for GCD over Q_{D}[X]. In this lecture, special techniques will be described for three problems:

(i) Computing determinants. The method to be described has similarities to the subresultant PRS algorithm, and in the modern form, is due to Bareiss [2, 3]. Bareiss noted that the method is known to Jordan. This method seems quite effective when D is the integers or univariate polynomials [5] (see also Sasaki and Murao [15]).

(ii)&(iii) Computing Hermite and Smith normal forms of integer matrices. These have applications to lattice-theoretic questions, solving linear Diophantine equations, and finitely generated Abelian groups.

The results for computing determinants apply to any domain D. For the Hermite and Smith normal forms, we describe the results for D = ℤ although the basic method could be extended to any UFD.

The set of m × n matrices with entries in D is denoted D^{m×n}. The (i, j)th entry of a matrix M is denoted (M)_{i,j} or (M)_{ij}.

§1. Sylvester’s Identity

Bareiss’ algorithm for computing determinants of matrices over D is based on a determinantal identity which we derive in this section. When applied to D = ℤ, Bareiss’ algorithm is polynomialtime.

It is instructive to understand why the usual Gaussian elimination is inadequate. Using fairly standard notations, suppose we initially have a matrix M^{(1)} with L-bit integers. In the kth stage, we transform M^{(k)} to M^{(k+1)}. The transformation applies to all entries of M^{(k)} with index (i, j) where i, j ≥ k. If the (i, j) entry of M^{(k)} is x_{ij}^{(k)}, we have

x_{ij}^{(k+1)} ← x_{ij}^{(k)} − x_{kj}^{(k)} \frac{x_{i1}^{(k)}}{x_{k1}^{(k)}} .

The entries in M^{(k+1)} are rational with bit size up to 4 times the bit sizes of entries in M^{(k)}. Thus after m steps, the entries can have bit size up to 4^{m}b. Pivoting does not seem to help (see Exercise). Hence new methods are needed. Edmonds [6] appears to be the first to give a polynomial time solution for this problem.

Let M ∈ D^{n×n} with (M)_{i,j} = x_{i,j}. Note that x_{ij} may also be a indeterminate if D is suitably defined. The (i, j)-cofactor of M is denoted [M]_{ij} and defined thus:

[M]_{ij} :=(−1)^{i+j} det M[i; j]  (1)

where M[i; j] is the matrix obtained by deleting the ith row and jth column of M. The adjoint of M is the matrix adj(M) whose (i, j)th entry is the (j, i)-cofactor of M (note the transposed subscripts). For any i, j = 1, …, n, it is easy to see that the sum

x_{i1}[M]_{j1} + x_{i2}[M]_{j2} + ⋯ + x_{in}[M]_{jn}

is equal to det M if i = j, and equal to 0 if i ≠ j. This immediately yields the following fundamental identity:

M · adj(M) = det(M) · I.  (2)

Taking determinants on both sides, it follows that det(adjM) = det(M)^{n−1}. If det(M) is an invertible element of D, we infer that the inverse of M exists and is given by

M^{−1} = (det(M))^{−1}adj(M).  (3)

Let ABC denote a triple matrix product of shape m × n × n × p (so B is an n-square matrix). For any b from domain D, we have the identity

bABC = A(bB)C.  (4)

This follows by looking at the typical element of bABC:

(bABC)_{rs} = b ∑_{i=1}^{n} a_{ri}(BC)_{is}

= b ∑_{i=1}^{n} a_{ri} ∑_{j=1}^{n}(B)_{ij}c_{js}

= ∑_{i=1}^{n} a_{ri} ∑_{j=1}^{n}(bB)_{ij}c_{js}

= ∑_{i=1}^{n} a_{ri}(bBC)_{is}
= (A(bB)C)_{rs}.

Next we express the matrix M in the form

M = [ \begin{matrix}A & B \\ C & D \end{matrix} ]

where A is (k − 1)-square and^{1} D is (n − k + 1)-square. For the following derivation, assume

δ = det A ≠ 0.

Lemma 1

δ^{n−k} det M = det(δD − C · adj(A)B).

Proof. If we express

M = [ \begin{matrix}A & 0 \\ C & I \end{matrix} ] [ \begin{matrix}I & A^{−1}B \\ 0 & D − CA^{−1}B \end{matrix} ]  (5)

then we see that

det M = δ det(D − CA^{−1}B)
δ^{n−k} det M = det(δD − δCA^{−1}B)
= det(δD − C · adj(A)B)

1Temporarily, “D” does not refer to the domain D.

where the last step exploits equations (3) and (4). Q.E.D.

But what is δD −C ·adj(A)B? To see this, introduce the “(r, s)-bordered matrix of order k” defined to be

M_{r,s}^{(k)} := [ \begin{matrix}x_{1,1} & x_{1,2} & ⋯ & x_{1,k−1} & x_{1,s} \\ x_{2,1} & x_{2,2} & ⋯ & x_{2,k−1} & x_{2,s} \\ ⋯ & & ⋯ & & ⋯ \\ x_{k−1,1} & x_{k−1,2} & ⋯ & x_{k−1,k−1} & x_{k−1,s} \\ x_{r,1} & x_{r,2} & ⋯ & x_{r,k−1} & x_{r,s} \end{matrix} ]

for k ≤ min{r, s}. By definition, M_{r,s}^{(1)} is the 1 × 1 matrix [x_{r,s}]. Also, define

x_{rs}^{(k)} := det M_{rs}^{(k)}.  (6)

For instance, M = M_{nn}^{(n)} and det M = x_{nn}^{(n)}. Now let us look at a typical element of δD−C ·adj(A)B: for r ≥ k and s ≥ k, we have

(δD − C · adj(A)B)_{r−k+1,s−k+1} = δx_{r,s} − (C · adj(A)B)_{r−k+1,s−k+1}
= δx_{r,s} − ∑_{i=1}^{k−1} C_{r−k+1,i} ∑_{j=1}^{k−1}(adj(A))_{ij}B_{j,s−k+1}

= δx_{r,s} − ∑_{i=1}^{k−1} x_{r,i} ∑_{j=1}^{k−1}[A]_{ji}x_{j,s}

where [A]_{ji} is the (j, i)-cofactor of A. But the last expression can be seen to be equal to the determinant of the bordered matrix M_{rs}^{(k)} (cf. exercise below giving the cofactors of x_{r,i}x_{j,s}). This proves:

Lemma 2

(δD − C · adj(A)B)_{rs} = x_{rs}^{(k)}.

Note that δ = x_{k−1,k−1}^{(k−1)}. Combining the last two lemmas:

Lemma 3 (Sylvester’s identity)

(x_{k−1,k−1}^{(k−1)})^{n−k} det M = det [ \begin{matrix}x_{k,k}^{(k)} & x_{k,k+1}^{(k)} & ⋯ & x_{k,n}^{(k)} \\ x_{k+1,k}^{(k)} & x_{k+1,k+1}^{(k)} & ⋯ & x_{k+1,n}^{(k)} \\ ⋮ & & ... & ⋮ \\ x_{n,k}^{(k)} & x_{n,k+1}^{(k)} & ⋯ & x_{n,n}^{(k)} \end{matrix} ] .

Exercises

Exercise 1.1: (i) (Wilkinson 1961) With the notations for Gaussian elimination in the introduction, let us now assume total pivoting. Show that

|a_{ij}^{(k)}| ≤ k^{1/2} (2 · 3^{1/2} · 4^{1/3} · · · · · k^{1/(k−1)})^{1/2} ∥A∥_{∞}.

§2. Determinant Computation  Lecture X  Page 4

Note: a rough estimate is that log (2 · 3^{1/2} · 4^{1/3} · · · · · k^{1/(k−1)}) is O(log^{2} k). Thus the magnitude of the entries are polynomial. (ii) Why is it not obvious that this leads to a polynomial time solution for integer matrices? (iii) (Open) Construct examples with exponential bit sizes in the intermediate entries. (You may assume that no pivoting is needed.) □

Exercise 1.2: (Bareiss) Show that Sylvester’s identity holds even when x_{k−1,k−1}^{(k−1)} = 0. HINT: perturb the singular submatrix. □

§2. Fraction-free Determinant Computation

Lemma 3 with n − k = 1 is called the “first order” Sylvester identity. With the notations of the previous section, this identity for the matrix M_{i,j}^{(k+1)} amounts to

x_{k−1,k−1}^{(k−1)} det M_{ij}^{(k+1)} = det [ \begin{matrix}x_{k,k}^{(k)} & x_{k,j}^{(k)} \\ x_{i,k}^{(k)} & x_{i,j}^{(k)} \end{matrix} ] .

Hence,

x_{ij}^{(k+1)} = \frac{x_{k,k}^{(k)}x_{i,j}^{(k)} − x_{i,k}^{(k)}x_{k,j}^{(k)}}{x_{k−1,k−1}^{(k−1)}} .  (7)

The important point is that the division by x_{k−1,k−1}^{(k−1)} in this equation is exact (i.e., with no remainder). Equation (7) is the defining step of the fraction-free Gaussian elimination algorithm of Bareiss (and Jordan):

Bareiss Algorithm

Input: M an n-square matrix,

assuming its principal minors x_{kk}^{(k)} are all non-zero.

Output: The matrix entry (M)_{n,n} contains the determinant of M. In general, for i ≥ k, we have (M)_{ik} = x_{ik}^{(k)}, (M)_{ki} = x_{ki}^{(k)}.

1. (M)_{0,0} ← 1; {Note: (M)_{0,0} is a special variable.}
2. for k = 1, …, n − 1 do
3. for i = k + 1, …, n do
4. for j = k + 1, …, n do
5. (M)_{ij} ← \frac{(M)_{ij}(M)_{kk}−(M)_{ik}(M)_{kj}}{(M)_{k−1,k−1}}

The program structure of this algorithm amounts to a simple triple-loop, as in the standard Gaussian elimination. Its correctness is easily shown by induction on k, and by appeal to equation (7) (we leave this as an exercise).

In case the assumption about principal minors turns out to be false, this is easily detected and the algorithm may be aborted. Alternatively, it is not hard to add the code (between lines 2 and 3) to perform some kind of pivoting: say, if (M)_{k−1,k−1} = 0 and some (M)_{k−1,i} ≠ 0 (i = k, …, n) then we can exchange the ith column with the k − 1st column. But we defer a more complete discussion of this to an extension of Bareiss’ algorithm below.

The division in line 5 is exact since (M)_{k−1,k−1} = x_{k−1,k−1}^{(k−1)}. Hence all computed values remain inside the domain D.

This is an “in-place” algorithm that destroys the contents of the original matrix. But we can easily preserve the original matrix is desired. The output M has the following shape:

M = [ \begin{matrix}x_{1,1}^{(1)} & ⋯ & & & & x_{1,n}^{(1)} \\ ⋮ & ... \\ & & x_{kk}^{(k)} & x_{k,k+1}^{(k)} & ⋯ & x_{k,n}^{(k)} \\ & & x_{k+1,k}^{(k)} & ... \\ & & ⋮ & ... \\ x_{n,1}^{(1)} & & x_{n,k}^{(k)} & & & x_{n,n}^{(n)} \end{matrix} ]

In other words, for each k = 1, …, n, there is an (rotated) L-shaped band in M that contains determinants of order k bordered matrices, as indicated.

In view of the definition of x_{ij}^{(k)} (equation (6)) as subdeterminants of M, we obtain at once:

Lemma 4 The intermediate values encountered in the algorithm has absolute values at most n^{n}2^{Ln} where 2^{L} bounds the absolute values of the entries of M.

Since the algorithm takes O(n^{3}) arithmetic steps, and each entry has at most n(log n + L) bits, we conclude that the bit-complexity of this method is

O(n^{3}M_{B}(n(log n + L)))

where M_{B}(s) is the bit-complexity of multiplying two s-bit integers.

One can exploit the higher order Sylvester identities to obtain analogous algorithms. We will not explore this but see Bareiss [2] for an analysis of an algorithm exploiting the second order identity. We will instead describe two other extensions.

Extended Bareiss Algorithm. We extend the algorithm to matrices of arbitrary m × n shape and of general rank ρ. This can be formulated as follows:

(*) Given a matrix A ∈ ℤ^{m×n}, we want to compute its rank ρ, a permutation P of its rows, a permutation Q of its columns, and non-zero values d_{1}, …, d_{ρ}, such that each d_{i} is the ith principal minor of P AQ.

This will be needed for the Hermite normal form algorithm in §7.

Partially Processed Matrices. To describe the solution, it is useful to introduce some terminology. Let A, M ∈ ℤ^{m×n}. We say that M is a partially processed version of A if each (i, j)th entry (M)_{i,j} is an (i, j)-bordered determinant of A of some order k = k(i, j). Clearly we have

1 ≤ k ≤ min{i, j}.

If k(i, j) = min{i, j} (respectively, k(i, j) = 1), we say the (i, j)th entry is fixed (resp., original). We call k(i, j) the fixing order of the (i, j)th entry. For instance, A is a partially processed version of itself, with every entry original. If every entry of M is fixed (for A), then we say M is said to be completely processed (for A). In this terminology, we can view Bareiss’ algorithm on input A as trying to fix every entry of A. An operation of the form

(M)_{i,j} ← \frac{(M)_{k,k}(M)_{i,j} − (M)_{i,k}(M)_{k,j}}{(M)_{k−1,k−1}}  (8)

a fixing step for the (i, j)th entry. The fixing step is valid provided the fixing order of (M)_{k,k}, (M)_{i,m}, (M)_{i,k} and (M)_{k,m} is k and the fixing order of (M)_{k−1,k−1} is k − 1. If M is partially processed before a valid step for the (i, j)th entry, it remains partially processed after the step, with the fixing order of the (i, j) entry equal to k + 1 (cf. equation (10))

Let us return to problem (*). Suppose M initialized to A and the row and column permutations P, Q are initially identities. Inductively, assume M is a partially processed version of P AQ. We proceed in stages. In stage s (s = 1, …, ρ) our goal is ensure the following properties:

(i) Every entry in the sth principal subdeterminant of M is fixed. (ii) The first s diagonal entries of M are non-zero. (iii) To maintain an index m_{0} ≥ s such that each row i > m_{0} is known to be dependent on the first s − 1 rows. Moreover, the entries in rows s + 1 to m_{0} are original. Initially, m_{0} = m.

For s ≥ 1, suppose that stage s − 1 is completed. [For s = 1, the properties (i)-(iii) are vacuously true.] Observe that we can fix the first s − 1 entries of the sth row of C in O(s^{2}) fixing steps. Applying this observation once more, we can fix the first s entries in the sth column in O(s^{2}) fixing steps. Thus goal (i) is attained. The problem is that the (s, s)th entry may be zero. We may go on to fix the first s entries of column j for j = s + 1, s + 2, …, n. There are two cases: (A) We find a column j which can be exchanged with column s to satisfy goal (ii). Then we exchange column j and column s, update permutation Q. If s = m_{0}, we halt, else go to stage s + 1. (B) No such column exists. In this case we conclude that row s is dependent on the first s − 1 rows. If m_{0} = s, we stop, since ρ = s − 1. Otherwise, we exchange rows m_{0} with row s, decrement m_{0}, update permutation P and continue with stage s.

This completes the description of our solution to problem (*). For reference, call this the Extended Bareiss Algorithm. It’s correctness is clear. For its complexity, first note that each matrix entry is a minor of A and hence has bit-size at most L′:= ρ(lg ρ + L) where L = lg ∥A∥_{∞}. Second, we use O(mnρ) arithmetic operations on integers of bit-size O(ρ(lg ρ + L)). This is because the fixing step for each entry is applied at most ρ + 1 times. To see this, we verify two facts: (a) In case (A), when we exchange columns s and j, note that the first s entries in the columns between s and j have been fixed. This information is still valid after the column exchange. So we just need a Boolean flag at each entry to indicate whether it is fixed or not. (b) In case (B), note that after the exchange between rows s and row m_{0}, entire row m_{0} is already fixed and row s is original. Thus we have:

Theorem 5 On input A ∈ ℤ^{m×n} with rank ρ ≥ 1, the Extended Bareiss Algorithm has bit complexity O(mnρM_{B}(L′)) where L′ = ρ(lg ρ + L) and L = lg ∥A∥.

A Preprocessing Problem. We note another useful extension to Bareiss’ algorithm. If N ∈ ℤ^{m×(m−1)}, let N_{i} (i = 1, …, m) be the submatrix of N with the ith row deleted. Consider this problem:

(**) Given N ∈ ℤ^{m×(m−1)}, compute the m subdeterminants

det N_{1}, …, det N_{m}.  (9)

This is a “preprocessing problem” in the sense that once the det N_{i}’s are available, for any given column a = (a_{1}, …, a_{m})^{T} ∈ ℤ^{m}, if A is obtained by appending a to N, we can quickly compute det A as

det A = ∑_{i=1}^{m} a_{i}(−1)^{i+m+1} det N_{i}.

The solution is simple: let M be the m × m matrix obtained from N by appending a column

Z = (Z_{1}, …, Z_{m})^{T}

where the Z_{i}’s are indeterminates. Let X = (x_{i,j})_{i,j=1}^{m,m} be the output matrix when M is input to Bareiss’ algorithm. The entries x_{i,j} are computed as expected. But the entry x_{i,m} in the last column is a linear function in Z_{1}, …, Z_{i}, which we denote by f_{i}(Z_{1}, …, Z_{i}).

In the notations of §1, let x_{i,j}^{(k)} denote the (i, j)-bordered determinant of M of order k (i, j ≥ k). Thus each output entry x_{i,j} is equal x_{i,j}^{(k)} where k = min{i, j}. Let δ_{k} = x_{k,k}^{(k)}. For example, the reader may verify that

f_{1} = Z_{1},
f_{2} = δ_{1}Z_{2} − x_{2,1}Z_{1},
f_{3} = δ_{2}Z_{3} − x_{3,2}Z_{2} − (\frac{δ_{2}x_{3,1} − x_{3,2}x_{2,1}}{δ_{1}} ) Z_{1}.

For any (a_{1}, …, a_{m}) ∈ ℤ^{m}, the value of f_{m}(a_{1}, …, a_{m}) yields the determinant of the matrix M where each Z_{i} is replaced by a_{i}. Thus, up to signs, the desired subdeterminants (9) may be read off from the coefficients of f_{m}.

What is the complexity of this procedure? The bit-sizes of entries in the first m − 1 columns of M and the time to compute them are exactly as in Bareiss’ algorithm. Consider the mth column. In stage k + 1 (k = 1, …, m − 1) of the outermost for-loop in Bareiss’ algorithm, the entries of the column m that are computed are x_{i,m}^{(k+1)} (i = k + 1, …, m). We have

x_{i,m}^{(k+1)} = \frac{x_{k,k}^{(k)}x_{i,m}^{(k)} − x_{i,k}^{(k)}x_{k,m}^{(k)}}{x_{k−1,k−1}^{(k−1)}}  (10)

= \frac{δ_{k}x_{i,m}^{(k)}(Z_{1}, …, Z_{k−1}, Z_{i}) − x_{i,k}^{(k)}f_{k}(Z_{1}, …, Z_{k})}{δ_{k−1}}.  (11)

The linear functions x_{i,m}^{(k)}(Z_{1}, …, Z_{k−1}, Z_{i}) and f_{k}(Z_{1}, …, Z_{k}) each has k coefficients that are minors of N of order k −1, and hence has bit-size at most m(L+lg m). Hence it takes takes O(k) arithmetic operations to compute x_{i,m}^{(k+1)}. Summing over the cost for computing the entries of column m, we again have O(m^{3}) arithmetic operations on integers of bit size O(m(L + lg m)). So the overall complexity is exactly as in the original Bareiss’ algorithm.

Exact Division. Exact division turns out to be slightly more efficient than division-with-remainder (by a small constant factor). We briefly describe the method (see Jebelean [11] for more details). Suppose C = A·B is an integer product. Consider the problem of computing B given C, A where integers are in base β. Let 0 ≤ ℓ(A) < β denote the least significant digit (LSD) of A. Then it is easy to see that ℓ(C) ≡ ℓ(A)ℓ(B)(mod β). Thus

ℓ(B) = (ℓ(C) · (ℓ(A)^{−1}) mod β)

provided ℓ(A) is invertible mod β. For simplicity, assume β is prime so that this is always possible. This immediately leads to an exact division algorithm that produces the digits of B sequentially, beginning with the LSD of B. Clearly, this is the opposite of the classical division algorithm [12], and avoids the “guess-and-correct” step of the classical method:

Exact Division Input: A, C ∈ ℤ in a prime base β, A|C. Output: B such that AB = C.

NORMALIZATION:

1. while β|A do

2. A ← A/β; C ← C/β.

MAIN LOOP:

3. while C > 0 do

4. b ← ℓ(C)/ℓ(A) mod β;

5. Output b;

6. C ← (C − b · A)/β.

Let len(A) denote the number of digits in A. Step 6 is considered the inner loop. To speed up this step, observe that only the lowest len(B) digits of C are involved in the inner loop. Hence the main loop can be improved as follows:

⋯

MAIN LOOP:

3. for k ← (len(C) − len(A) + 1) downto 1 4. b ← ℓ(C)/ℓ(A) mod β;

5. Output b;

6. C ← ((C − b · A) mod β^{k})/β.

If β is a power of 2, then ℓ(A) would be invertible if A is odd. We achieve this by a simple modification to the normalization stage, namely, by inserting steps 2.1 and 2.2 below:

NORMALIZATION:

1. while β|A do

2. A ← A/β; C ← C/β. 2.1 while 2|A do 2.2 A ← A/2; C ← C/2. MAIN LOOP:

⋯

At the end of normalization, A is odd, and hence ℓ(A) is odd. This guarantees that ℓ(A) is invertible. The bit analysis of this algorithm is left to an exercise.

§3. Matrix Inversion  Lecture X  Page 9

Exercises Exercise 2.1: M is a square matrix with (M)_{ij} = x_{ij} for all i, j. The “cofactor” of x_{ij}x_{i′j′} is

defined to be the expression E that is multiplied by x_{ij}x_{i′j′} when we collect terms in the determinant of M that involve both x_{ij} and x_{i′j′}. E.g., if M is 2 × 2, the cofactor of x_{12}x_{21} is −1 and the cofactor of x_{11}x_{22} is 1. If i = i′ or j = j′ then E = 0; otherwise, show that

E = (−1)^{i+i′+j+j′+δ} det M[i, i′; j, j′]

where

δ = \{ \begin{matrix}1 & if (i > i′) ⊕ (j > j′), (⊕ is exclusive-or) \\ 0 & else, \end{matrix}

and M[i, i′; j, j′] is the submatrix of M obtained by deleting rows i and i′ and deleting columns j and j′. □

Exercise 2.2: Show that adj(adjA) = det(A)^{n}A. □

Exercise 2.3: What is the 3 × 3 matrix analogue of equation (5)? □

Exercise 2.4: Carry out the d-th order version of Bareiss algorithm, by exploiting the order d Sylvester identity. For instance, for d = 2, we must construct x_{ij}^{(k)} for even values of k, evaluating 3 × 3 determinants. □

Exercise 2.5: Modify Bareiss’ algorithm for efficient computation of a matrix with rational entries. Carry out comparison experiments in this setting of rational inputs. □

Exercise 2.6: Suppose M is n × m where n ≤ m. In Bareiss’ algorithm, we replace line 4 with “for j = k +1, …, m do” (i.e., we extend change the limit from n to m). Describe the contents of the entries (M)_{n,j} for j = n, n + 1, …, m. What is the complexity of this modification? How is this related to the determinantal polynomials (Lecture III)? □

Exercise 2.7: In the exact division algorithm, show that when the length m of C is less than twice the length n of A, this method uses no more than half the number of bit-operations of the classical method. Quantify this bit-operation count more generally in terms of m, n. □

§3. Matrix Inversion

Matrix inverse is easily computed using the standard Gaussian elimination procedure. Following [7], let us compute the more general product

CA^{−1}B

where the triple product CAB has shape m × n × n × p. We proceed as follows: apply Gaussian elimination to the (m + n) × (n + p) matrix

M = [ \begin{matrix}A & B \\ −C & 0 \end{matrix} ] ,  (A is nonsingular)  (12)

and obtain (after operations that zero the first n columns below A)

M′ = [ \begin{matrix}A′ & B′ \\ 0 & D′ \end{matrix} ] .  (13)

Here 0 denotes a matrix of zeros. Note that if A is singular, we would have discovered this fact during Gaussian elimination. Henceforth, assume A is non-singular. We claim that D′ is precisely what we wanted:

Lemma 6

D′ = CA^{−1}B.

Block Gaussian elimination. This lemma is a slight generalization of lemma 1 to the non-square case. It is instructive to briefly consider Gaussian elimination for block-size elements. Let

M = [ \begin{matrix}A_{11} & A_{12} & ⋯ & A_{1n} \\ A_{21} & A_{22} & ⋯ & A_{2n} \\ ⋮ & & ... & ⋮ \\ A_{m1} & A_{m2} & ⋯ & A_{mn} \end{matrix} ]

where A_{ij} is a k_{i} × ℓ_{j} matrix (“block”). We may say that M has shape (k_{1}, …, k_{m}) × (ℓ_{1}, …, ℓ_{n}). Consider the following transformation of M: for 1 ≤ i, j ≤ m, i ≠ j, and any ℓ_{i} × ℓ_{j} matrix H, we replace the ith row of M with the ith row plus H left multiplied by the jth row. That is,

A_{ik} ← A_{ik} + HA_{jk}, (for k = 1, …, n).

This is equivalent to left multiplying M by the matrix

T_{i,j}(H) := [ \begin{matrix}I_{1} & ⋯ & 0 & ⋯ & 0 & ⋯ & 0 \\ ⋮ & ... & ⋮ & & ⋮ & & ⋮ \\ 0 & ⋯ & I_{i} & ⋯ & H & ⋯ & 0 \\ ⋮ & & ⋮ & ... & ⋮ & & ⋮ \\ 0 & ⋯ & 0 & ⋯ & I_{j} & ⋯ & 0 \\ ⋮ & & ⋮ & & ⋮ & ... & ⋮ \\ 0 & ⋯ & 0 & ⋯ & 0 & ⋯ & I_{m} \end{matrix} ]

where I_{j} is the ℓ_{j}-square identity matrix, the I_{j}’s occur along the main diagonal, 0-blocks occur everywhere else except at the (i, j)th position which is occupied by H. Clearly, for any minor Δ of M (now viewed as an ordinary matrix of shape (∑_{i=1}^{m} k_{i}) × (∑_{j=1}^{n} ℓ_{j})) that contains rows in the ith and jth block rows of M, this operation preserves Δ. If m = n and for all i, k_{i} = ℓ_{i} and A_{i,i} is non-singular, then we can effect Gaussian elimination at this block level: for instance, T_{i1}(A_{i1}A_{11}^{−1}) will make the (i, 1)th entry zero. In particular, if A is square and non-singular, we may transform the matrix

N = [ \begin{matrix}A & B \\ C & D \end{matrix} ]  (14)

to

N′ = [ \begin{matrix}A′ & B′ \\ 0 & D − CA^{−1}B \end{matrix} ] .

We may conclude:

Lemma 7 The rank of N is equal to the rank of A iff D = CA^{−1}B.

We complete the proof of lemma 6. Suppose, instead of M in equation (12), we had

P = [ \begin{matrix}A & B \\ −C & −CA^{−1}B \end{matrix} ] .

§4. Hermite Normal Form  Lecture X  Page 11

Then by lemma 7, the rank of P is equal to the rank of M. Applying Gaussian elimination to P to zero the entries below A, we obtain

P′ = [ \begin{matrix}A′ & B′ \\ 0 & D′ − CA^{−1}B \end{matrix} ] .

The matrices A′, B′, D′ here are the same as those in M′ (equation (13)) since the same multipliers were used to do the elimination. But the rank of P′ equals the rank of P and hence of M. This can only mean that D′ − CA^{−1}B = 0, as we wanted shown.

Exercises

Exercise 3.1: Let N be square in equation (14).

(i) Show that det N = det(AD − ACA^{−1}B) if A^{−1} exists.

(ii) Under what conditions does det N = det(AD − CB)? det N = det(AD − BC)?

(iii) Consider the four cases of these identities depending on which of the 4 blocks of N are non-singular. (iv) Obtain the corresponding formulas using column operations (multiplication from the right). □

Exercise 3.2: Use the method to compute CA^{−1}B to solve the system Ax = b of linear equations. □

Exercise 3.3: For n > 1, let U be the n-square matrix all of whose entries are 1 and let S = U − I. So S has zeros along its main diagonal. Compute S^{−1} using the above algorithm for small n. Note: S^{−1} = \frac{1}{n−1}U − I. □

Exercise 3.4: Let T (n) be the algebraic complexity of computing the determinant of an n-square matrix. Show that T (n) = O(MM(n)) where MM(n) is the complexity of multiplying two n-square matrices. □

§4. Hermite Normal Form

Let A, B ∈ ℤ^{m×n}. In §VIII.1, we viewed the columns of A as a generating set of the lattice Λ(A) ⊆ ℤ^{m}. Alternatively, we may regard Λ(A) as a subgroup of the Abelian group ℤ^{m}. We have shown (§VIII.1)

Λ(A) = Λ(B) iff A = B · U  (15)

for some integer unimodular matrix U. The original proof requires A and B to be bases but this assumption can be dropped (see below). This raises the question: how can we decide if two matrices A, B generates the same subgroup (or lattice)? The result (15) does not appear helpful for this purpose – there is no obvious way to find U even if one is known to exist. In this lecture, and unlike §VIII.1, we will no longer assume that A have rank n, so the columns of A forms a generating set but not necessarily a basis of Λ(A). The tool for answering such questions about Λ(A) is a normal form which we now describe.

By applying the elementary integer column operations (§VIII.1) to A, we can transform A to a matrix H = H(A) of the following shape:

1. For some r = 1, …, n, the first r columns of H are non-zero and the remaining columns are zero.

2. For i = 1, …, r, let (H)_{c(i),i} be the first non-zero entry of column i. Then

1 ≤ c(1) < c(2) < ⋯ < c(r) ≤ n.  (16)

3. (H)_{c(i),i} > 0.

4. If 1 ≤ j < i ≤ r then 0 ≤ (H)_{c(i),j} < (H)_{c(i),i}.

Such a matrix H = H(A) is said to be in Hermite normal norm (HNF), or H is the HNF of A. For instance, the following 6 × 4 matrix H_{1} is in HNF:

H_{1} = [ \begin{matrix}2 & 0 & 0 & 0 \\ −1 & 0 & 0 & 0 \\ 0 & 3 & 0 & 0 \\ 0 & 1 & 2 & 0 \\ 2 & 10 & −4 & 0 \\ 9 & −3 & 0 & 0 \end{matrix} ] .  (17)

If only the first two conditions in the definition of HNF hold, H is said to be in column echelon form. Column echelon form is just a generalization of “lower triangular matrices”. If H is a column echelon form of A, then rows c(1), …, c(r) are called the critical rows of H. Each of the entries (H)_{c(i),i} is called a critical entry of H. The number r of non-zero columns in H is called the rank of A. For the matrix H_{1} in (17), the critical rows are the 1st, 3rd and 4th, the critical entries are (1, 1), (3, 2) and (4, 3).

Since the elementary column operations preserve the underlying lattice, Λ(A) = Λ(H(A)). It is then easy to see: row i is critical in H(A) iff there is a vector ξ ∈ Λ(H(A)) whose first non-zero entry is in the ith position. This latter characterization of critical rows depends only on the lattice. This shows that the set of critical rows of H(A) depends only on Λ(A).

Theorem 8 Let A, B have the same shape. Λ(A) = Λ(B) if and only if A and B has the same Hermite normal form: H(A) = H(B).

Proof. If H(A) = H(B) then clearly Λ(A) = Λ(B). In the other direction, let us assume Λ(A) = Λ(B). We wish to prove H(A) = H(B). First observe that A and B has the same rank, d. This is because A, B have the same set of critical rows, by a preceding observation. But d is the number of critical rows.

We use induction on d. We assume that d ≥ 1 since d = 0 is trivial. Let a, b be the first columns of H(A) and H(B), respectively. Let c(1) and c′(1) denote the first non-zero entry of a and b, respectively. It is easy to see that c(1) = c′(1). Indeed, these two entries must be identical since they each generate all the c(1)th entries of vectors in Λ(A) = Λ(B).

Let H_{1}(A) and H_{1}(B) be obtained by deleting a and b, respectively. Notice that Λ(H_{1}(A)) generates the subgroup of Λ(A) whose first c(1) entries are zero. Similarly Λ(H_{1}(B)) generates the subgroup of Λ(B) whose first c(1) entries are zero. Since Λ(A) = Λ(B), we conclude that Λ(H_{1}(A)) = Λ(H_{1}(B)). Since H_{1}(A) and H_{1}(B) are in HNF, we deduce by induction that H_{1}(A) = H_{1}(B).

It remains to prove that a = b. Suppose a − b is not identically zero. It follows from the preceding that the first c(1) entries of a−b are zero. If k is the first non-zero column of a−b, then k > c(1) and it follows that there exists a column c in H_{1}(A) whose first non-zero entry is in position k. If (a)_{k} denotes the kth component of a, then by the definition of HNF, (c)_{k} > (a)_{k} ≥ 0 and (c)_{k} > (b)_{k} ≥ 0. Hence (a − b)_{k} has absolute value less than (c)_{k}. This is impossible since a − b ∈ Λ(H_{1}(A)) means that (a − b)_{k} must be a multiple of (c)_{k}. Q.E.D.

As corollary, we also see that (15) holds without the restriction that the columns of A, B (respectively) are linearly independent. It suffices to show that if Λ(A) = Λ(B) then there is a unimodular matrix U such that AU = B. But there are unimodular matrices U_{A}, U_{B} such that AU_{A} = H(A) = H(B) = BU_{B}. Take U = U_{A}U_{B}^{−1}.

It follows that our basic question of deciding if Λ(A) = Λ(B) is reduced to computing and comparing their HNF’s. Other questions such as checking whether a given vector ξ belongs to Λ(A) can similarly be answered (Exercise). We next address the computation of HNF.

Generic HNF algorithm. Assume the input is the m × n matrix A. By the “subrow” at an entry (A)_{i,j} we mean the set of entries of the form (A)_{i,k} for k = j, j + 1, …, n. The best way to understand this algorithm is to imagine that our task (the main loop below) is to determine the critical rows of A.

Generic HNF Algorithm:

Input: an m × n integer matrix A.

Output: the Hermite normal norm of A. MAIN LOOP:

1. Initialize i ← 1 and j ← 1.

2. while i ≤ m and j ≤ n do:

3. While the subrow at (A)_{ij} has only zero entries, increment i.

4. Now assume the subrow at (A)_{ij} has non-zero entries.

5. By adding multiples of one column to another, eliminate all but one non-zero element in the subrow at (A)_{ij}.

6. By a column-exchange, bring this sole non-zero entry to position (i, j) and increment j.

end{while}

CLEAN UP:

7. At this point, the matrix is in column-echelon form.

8. By adding multiples of one column to another, achieve the remaining two conditions for the definition of HNF.

Exercise 4.1: (i) There are several ways to fill in details in the generic algorithm. Describe a reasonable choice. (ii) For your version of the generic algorithm, analyze the number of arithmetic operations for a m × n input matrix whose entries are L-bit integers. (iii) Bound the bit complexity of your algorithm. □

Although part (i) of the exercise should give a polynomial bound in m, n, L, the bit complexity in part (ii) should be exponential in min{m, n}. Note that these exponential bounds arise from potential sizes of intermediate matrix entries; the final entries in the HNF can be shown to be polynomially bounded (below). Nevertheless, it is an open problem to construct examples that actually exhibit exponential behavior. In random examples, huge intermediate entries routinely appear. For instance, Hafner and McCurley [8] reported that for random 20 × 20 matrices with entries between 0 and 10, most gave rise to an entry exceeding 10^{500}. One example has an entry exceeding 10^{5011}. We will develop a modular technique to achieve polynomial bit-complexity bounds.

Invariants of the Hermite normal form. Let A ∈ ℤ^{m×n} and 1 ≤ k ≤ min{m, n}. For 1 ≤ i_{1} < i_{2} < ⋯ < i_{k} ≤ m, let A(i_{1}, i_{2}, …, i_{k}) denote the submatrix of A formed by rows i_{1}, …, i_{k}. Let γ(A; i_{1}, …, i_{k}) denote the GCD of all the k × k subdeterminants (i.e., order k minors) of A(i_{1}, i_{2}, …, i_{k}). Note that γ(A; i_{1}, …, i_{k}) = 0 iff every order k minor of A(i_{1}, i_{2}, …, i_{k}) is zero.

Lemma 9 The elementary integer column operations on A preserve γ(A; i_{1}, i_{2}, …, i_{k}).

Proof. The column operations that interchange two columns or multiply a column by −1 does not change the GCD (since GCD is defined up to associates and we always pick the positive member.) Suppose c ∈ ℤ times the ith column of A is added to the jth column. Certain of the k × k subdeterminants of A(i_{1}, i_{2}, …, i_{k}) changes. If a subdeterminant value D is changed, say to D′, it is easy to see that D − D′ = ±c · D′′ where D′′ is another subdeterminant. Moreover, D′′ is among the subdeterminants of A(i_{1}, …, i_{k}) that has not changed. To see this, observe that a subdeterminant D is changed iff D involves column j but not column i. Schematically, if the old GCD is GCD(…, D, …, D′′, …) then the new one is

GCD(…, D′, …, D′′, …) = GCD(…, D ± c · D′′, …, D′′, …).

But the later expression is equal to the old GCD. Q.E.D.

Let r be the rank of A. For k = 1, …, r, we use the shorthand

γ_{k}(A) := γ(A; c(1), c(2), …, c(k)).

We define γ_{i}(A) = 0 for i = r + 1, …, n.

Corollary 10 The value γ_{k}(A) is invariant under the elementary column operations on A. If H is the Hermite normal form of A, then the product of the first k critical values of H is equal to γ_{k}(A). In particular, γ_{k}(A) divides γ_{k+1}(A) for k = 1, …, n − 1.

Exercises

Exercise 4.2: Discuss other strategies for implementing the generic HNF algorithm. □

Exercise 4.3: Solve the following problem efficiently:

(i) Given x ∈ ℤ^{m} and A ∈ ℤ^{m×n}, decide if x ∈ Λ(A). If x ∈ Λ(A), find the n-vector ξ such that Aξ = x. (ii) Given A, B ∈ ℤ^{m×n}, check whether Λ(A) = Λ(B). In case they are equal, construct matrices U, V such that A = BU and B = AV. [These matrices need not be unimodular unless A, B are bases for Λ(A).] □

Exercise 4.4: Suppose A, B are both in column-echelon form, and Λ(A) ⊆ Λ(B). If the critical entries in A and the critical entries in B are in the same position and corresponding critical entries have the same values, then Λ(A) = Λ(B) □

§5. GCD Bound and Algorithm  Lecture X  Page 15

Exercise 4.5: (i) Every subgroup of ℤ^{n} is finitely generated.

(ii) Every finitely generated Abelian group G is isomorphic to a subgroup of ℤ^{n} for some n. □

Exercise 4.6: Call H a generalized HNF for A if H is in HNF and obtained from A by the usual elementary column operations, but now we allow the permutation of rows as well. How do the various generalized HNF of A relate to each other? □

Exercise 4.7: (open) Given a matrix A ∈ ℤ^{m×n} and its Hermite normal form H(A). Let L bound the bit sizes of entries in A and H(A). What is best upper bound B(m, n, L) such that there exists a sequence of elementary integer column operations from A to H(A) where all intermediate entries have bit size at most B(m, n, L). □

§5. A Multiple GCD Bound and Algorithm

Computing multiple GCD’s is a key subproblem in Hermite normal forms. For example, in the generic HNF algorithm (§4), the process of zeroing out all but one entry of a subrow amounts to computing the multiple GCD of the entries in the subrow. Of course, multiple GCD can be reduced to simple GCD, i.e., GCD for two elements. But this is not the most efficient method. In this section, we present an efficient multiple GCD algorithm over integers. This is based on a co-factor bound that we first derive.

In the following, fix

(a_{1}, a_{2}, …, a_{k}) (k ≥ 2)

such that the a_{i}’s are distinct and positive. Let d = GCD(a_{1}, …, a_{k}). By definition, an integer sequence (s_{1}, …, s_{k}) is called a co-factor of (a_{1}, …, a_{k}) if

d = s_{1}a_{1} + s_{2}a_{2} + ⋯ + s_{k}a_{k}.

The “co-GCD problem” refers to the problem of computing a co-factor of (a_{1}, …, a_{k}). Note that once a co-factor is available, we can easily compute the GCD. Our first goal is to prove the existence of a co-factor with each |s_{i}| upper bounded by a_{1}. We use an argument of Hongwei Cheng:

Lemma 11 If d = GCD(a_{1}, …, a_{k}) then there exists a co-factor (s_{1}, …, s_{k}) for (a_{1}, …, a_{k}) such that

|s_{1}| ≤ \frac{a_{k}}{2} ,
|s_{i}| < \frac{a_{i−1} + a_{k}}{2}  i = 2, …, k − 1,
|s_{k}| ≤ \frac{d}{a_{k}} + \frac{a_{k−1}}{2} .

Proof. Suppose (t_{1}, …, t_{k}) is any co-factor for (a_{1}, …, a_{k}). Define

s_{i}:= \{ \begin{matrix}t_{i} − q_{i}a_{k} & i = 1, …, k − 1 \\ t_{k} + ∑_{j=1}^{k−1} q_{j}a_{j} & i = k, \end{matrix}

where q_{1}, …, q_{k} ∈ ℤ are to be specified. It is not hard to check that (s_{1}, …, s_{k}) is also a co-factor for (a_{1}, …, a_{k}). We now define the q_{i}’s inductively. Pick q_{1} to be the symmetric quotient (§II.3) of t_{1} divided by a_{k}. Then |s_{1}| ≤ a_{k}/2, as desired. Now consider the partial sums

S_{i} = ∑_{j=1}^{i} s_{j}a_{j}.

Thus S_{1} = a_{1}s_{1} and |S_{1}| ≤ a_{1}a_{k}/2. Inductively, assume S_{i−1} has been defined so that |S_{i−1}| ≤ a_{i−1}a_{k}/2. For i = 2, …, k − 1, define q_{i} so that |S_{i}| ≤ a_{i}a_{k}/2. This is clearly possible since S_{i} = S_{i−1} + a_{i}s_{i} = S_{i−1} + a_{i}(t_{i} − q_{i}a_{k}). It follows that

|a_{i}s_{i}| = |S_{i} − S_{i−1}| ≤ \frac{a_{i}a_{k}}{2} + \frac{a_{i−1}a_{k}}{2}
|s_{i}| < \frac{a_{k}}{2} + \frac{a_{i−1}}{2} ,

as desired. Finally, we bound |s_{k}|. By definition of S_{k}, we have S_{k} = d, the GCD of a_{1}, …, a_{k}. So |s_{k}a_{k}| = |S_{k} − S_{k−1}| ≤ d + \frac{a_{k−1}a_{k}}{2} or |s_{k}| ≤ \frac{d}{a_{k}} + \frac{a_{k−1}}{2} . Q.E.D.

Note that for k = 2, this lemma gives a well-known bound

|s_{1}| ≤ \frac{a_{2}}{2} , |s_{2}| ≤ 1 + \frac{a_{1}}{2} .

Suppose we further assume that

a_{1} > a_{2} > ⋯ > a_{k}.

Then we may further infer the bounds |s_{i}| ≤ a_{i−1} − (k − i + 1)/2 for i = 2, …, k and |s_{1}| ≤ a_{k}/2. This yields:

Corollary 12 For all a_{1} > a_{2} > ⋯ > a_{k} ≥ 2, there exists a co-factor (s_{1}, …, s_{k}) for (a_{1}, …, a_{k}) such that

∏_{i=1}^{k} |s_{i}| < \frac{a_{k}}{2} ∏_{i=1}^{k−1} (a_{i} − \frac{k − i + 1}{2}) ≤ ∏_{i=1}^{k}(a_{i} − 1).

This says that the output size of the co-GCD algorithm need not be larger than the input size.

We now address the question of computing a co-factor (s_{1}, …, s_{k}) satisfying the lemma. We will use a divide and conquer approach. We split a_{1}, …, a_{k} into two subsets according to the parity of their subscripts, and assume inductively that we have computed c, c′, t_{1}, …, t_{k} such that

c := GCD(a_{2}, a_{4}, …, a_{2⌊k/2⌋}) = ∑_{i=1}^{⌊k/2⌋} t_{2i}a_{2i},

c′ := GCD(a_{1}, a_{3}, …, a_{2⌊(k−1)/2⌋+1}) = ∑_{i=1}^{⌊(k−1)/2⌋} t_{2i+1}a_{2i+1},

By a call to the simple extended GCD on c, c′, we obtain d, t, t′ such that

d = tc + t′c′.

By induction, we may assume that |t_{i}| is bounded according to Lemma 11. In particular, |t_{i}| < a_{1} for i = 1, …, k. Similarly, |t| < a_{1} and |t′| < a_{1}. Define

s′_{1} = \{ \begin{matrix}t_{i}t & if i = even, \\ t_{i}t′ & if i = odd. \end{matrix}

Thus ∑_{i=1}^{k} s′_{i}a_{i} = d and |s′_{i}| < a_{1}^{2} for all i. Following the proof of Lemma 11, we may now reduce s′_{1}, …, s′_{k} to s_{1}, …, s_{k}:

Reduction Step:

1. s_{1} ← s′_{1} mod a_{k} and S_{1} ← s_{1}a_{1}.

2. for i = 2, …, k − 1 do

S_{i} ← (S_{i−1} + a_{i}s′_{i}) mod a_{i}a_{k};
s_{i} ← (S_{i} − S_{i−1})/a_{i};
q_{i} ← (s′_{i} − s_{i})/a_{k}.

3. s_{k} ← s′_{k} + ∑_{i=1}^{k−1} q_{i}a_{i}.

The mod operator here is the symmetric version (§II.3). Note that both divisions in step 2 are exact.

Let

L = lg max{a_{1}, …, a_{k}}.  (18)

The bit complexity of this reduction step is O(kM_{B}(L)). Let T (k, L) be the bit complexity of the overall recursive procedure. Clearly

T (k, L) = 2T (k/2, L) + O(kM_{B}(L) + M_{B}(L) log L)  (19)

where ‘kM_{B}(L)’ comes from the reduction step and ‘M_{B}(L) log L’ comes from the simple co-GCD computation for c, c′ (Lecture II). The solution is easily seen to be T (k, L) = O(k(log k + log L)M_{B}(L)). Thus we have:

Theorem 13 There is a multiple co-GCD algorithm with bit complexity

T (k, L) = O(k(log k + log L)M_{B}(L)).

On input (a_{1}, …, a_{k}), the output co-factor (s_{1}, …, s_{k}) satisfy the bounds of lemma 11.

Application to multiple LCM. There are many applications of multiple GCD. An obvious application is for computing the primitive factorization (§III.1) of an integer polynomial. Bareiss [2] states that “there is no question that for maximum efficiency in any integer-preserving elimination code, the elements of all the rows and columns respectively should be made relative prime to each other before starting the elimination process”. Here we consider another application: the computation of multiple LCM. Simple GCD and simple LCM are closely related: LCM(a, b) = ab/GCD(a, b). It is slightly more involved to relate multiple GCD with multiple LCM. Our multiple GCD algorithm computes intermediate information that can be used to obtain the multiple LCM. Specifically, on input (a_{1}, …, a_{k}) the information can be organized as a binary tree T of height ⌈lg k⌉ with k leaves, each associated with an a_{i}. At each internal node u of T, let S_{u} ⊆ {a_{1}, …, a_{k}} denote the a_{i}’s associated to the leaves of the subtree rooted at u. Assume we store at u the value GCD(S_{u}). It is then simple to extend the multiple GCD algorithm so that we recursively compute at node u the LCM of S_{u}. If v, w are the children of u, we use the formula

LCM(S_{u}) = LCM(LCM(S_{v}), LCM(S_{w})) = \frac{LCM(S_{v}) · LCM(S_{w})}{GCD(S_{u})}.

Let L′ = lg LCM(a_{1}, …, a_{k}). With L as in (18), we see that L′ ≤ kL. The cost of additional computation at u is O(M_{B}(L′)). Overall, the additional cost is O(kM_{B}(L′)). Combined with theorem 13, we conclude:

§6. Hermite Reduction Step  Lecture X  Page 18

Theorem 14 There is an algorithm to compute both the GCD and LCM of a_{1}, …, a_{k} with bit complexity O(k(log kM_{B}(L) + log L · M_{B}(L) + M_{B}(L′)) which is

O(k(log L · M_{B}(L) + kM_{B}(L)).

Remarks. Illiopoulous [10] obtains the co-factor bound |s_{i}| = O(a_{i}^{k}) by using a balanced binary tree and co-factor bounds for simple GCD’s. Lüneburg [13] gives the bound

|s_{i}| ≤ a_{1}/2

for all i except when i = i_{0} (for some i_{0}). Moreover, |s_{i_{0}}| ≤ (1 + a_{1}(k − 1))/2. The dependence on k in these bounds is somewhat unsatisfactory. Hongwei Cheng^{2} shows the uniform bound |s_{i}| < a_{1} for all i; lemma 11 follows his argument.

Exercises

Exercise 5.1: Suppose we apply lemma 11 with

a_{k−1} > a_{k−2} > ⋯ > a_{2} > a_{1} > a_{k}.

How does this compare to the bound in the corollary?

□

Exercise 5.2: Verify the solution to recurrence (19) by an appropriate induction. □

Exercise 5.3: Suppose a_{1}, …, a_{k} ∈ ℤ[X]. Give an efficient method for computing the extended GCD of a_{1}, …, a_{k}. □

Exercise 5.4: (Olydzko) Is there a function f(k) > 0 that goes to infinity as k → ∞ such that for all a_{1} > a_{2} > ⋯ > a_{k} ≥ 2, a cofactor (s_{1}, …, s_{k}) exists where |s_{i}| ≤ |a_{1}|/f(k)? □

§6. Hermite Reduction Step

Let a ∈ ℤ^{1×n} be a non-zero row n-vector, n ≥ 2. The goal in this section is to construct a unimodular matrix U ∈ ℤ^{n×n} such that a · U is zero except for its first entry. Since U is invertible, the non-zero entry of a · U must be equal to GCD(a). We will call the transformation

a ↦ U  (20)

a “Hermite reduction step problem”. The reason for this name is clear. For, if a is the first row in a m × n matrix, then A · U is essentially the first step of the generic HNF algorithm. Repeating this process suitably, we finally obtain a column echelon form of A, which is easily converted into the Hermite normal form.

First let us illustrate our approach for the case n = 2. Suppose a = (a_{1}, a_{2}) and we wish to find a 2 × 2 unimodular U such that aU = (g, 0) where g = GCD(a_{1}, a_{2}). So there exist integers s, t such that g = sa_{1} + ta_{2}. Moreover, s, t are relatively prime. Hence there exist integers s′, t′ such that ss′ + tt′ = 1. In fact, we may choose s′ = a_{1}/g and t′ = a_{2}/g. Then it is easy to see that

2Private communication, 1991.

(a_{1}, a_{2}) [ \begin{matrix}s & −t′ \\ t & s′ \end{matrix} ] = (g, g′)  (21)

for some g′. But g divides g′ (since g divides a_{1} and a_{2}). If e = g′/g then we see that the desired U can be taken to be P · P^{∗} where

P = [ \begin{matrix}s & −t′ \\ t & s′ \end{matrix} ] , P^{∗} = [ \begin{matrix}1 & −e \\ 0 & 1 \end{matrix} ] .

It turns out that U can always be written as the product of matrices P and P^{∗}. In fact, we will see that it is more useful to represent U as the pair (P, P^{∗}). We begin with a key lemma which shows how to construct P (see [9, p. 375] or [14]).

Lemma 15 Let u = (u_{1}, …, u_{n})^{T} be a column of integers with d = GCD(u). There exists an n × n integer matrix P with the following properties:

(i) The first column of P equals u.

(ii) det P = d.

(iii) ∥P ∥_{∞} < ∥u∥_{∞}.

Proof. We use induction on n. The result is trivial for n = 1 so let n ≥ 2. Let u′ = (u_{1}, …, u_{n−1})^{T} with d′ = GCD(u′). By induction, there is a matrix P′ with first column u′ and det P′ = d′. So there are integers s, t such that

sd′ + tu_{n} = d.  (22)

We may assume that |s| < |u_{n}| and |t| < |d′|. We claim (§5) that the desired matrix can be taken to be

P = [ \frac{u_{1}}{r} ; P′  \frac{u_{2}}{r} ; ⋮ ; \frac{u_{n−1}}{r} ; u_{n} 0 ⋯ 0 s ] ,

where r ∈ ℚ is to be determined. Part (i) is clearly satisfied. By expanding the determinant along the last row, we have

det P = s det P′ + u_{n} det P′′,

where P′′ is obtained by omitting the first column and last row of P. We want to choose r so that det P′′ = t (so that (ii) is satisfied). We observation that P′′ is rather similar to P′: the last column of P′′ is just 1/r times the first column of P′. In fact,

P′′ = P′ · C = P′ · [ \begin{matrix}0 & 0 & 0 & ⋯ & 0 & \frac{1}{r} \\ 1 & 0 & 0 & & 0 & 0 \\ 0 & 1 & 0 & & 0 & 0 \\ ⋮ & & & ... & & ⋮ \\ 0 & 0 & 0 & & 0 & 0 \\ 0 & 0 & 0 & ⋯ & 1 & 0 \end{matrix} ] .

Here C has a subdiagonal of 1’s and top-right entry of 1/r. Note that C simply circulates the columns of P′, moving the first column (multiplied by 1/r) to the last place and for the other columns, one place to the left. When n = 2, C is simply the matrix [\frac{1}{r}]. Then det P′′ = det C det P′ = (−1)^{n} \frac{1}{r}d′.

Thus part (ii) is satisfied with r = (−1)^{n}d′/t. To see that P is correct, note that the entries of P are integers since u_{i}/r = u_{i}t/d′ ∈ ℤ for i = 1, …, n − 1.

Finally, part (iii) claims that each entry of P is bounded by ∥u∥_{∞}. This is true of the entries of P′ (by induction) and also true of the non-zero entries in the rest of P, namely, s, u_{n} and u_{i}/r (in the last case, we use the fact that |t| < |d′|). Q.E.D.

Note that P is not unique since s, t are not unique. In our application, the vector u = (u_{1}, …, u_{n})^{T} satisfies GCD(u) = 1 so that the matrix P is unimodular.

Now we implement the above lemma.

Lemma 16 Let a = (a_{1}, …, a_{n}) ≠ 0 be an integer row vector and log ∥a∥_{∞} ≤ L.

(i) We can construct the unimodular matrix P of the previous lemma 15 in bit complexity

O(n(n + log L)M_{B}(L)).  (23)

(ii) Suppose b = (b_{1}, …, b_{n}) ∈ ℤ^{n} satisfies log ∥b∥_{∞} ≤ L then we can compute b · P from b, P in bit

complexity

O(nM_{B}(L + log n)).  (24)

Proof. (i) Let u = (u_{1}, …, u_{n})^{T} such that ∑_{i=1}^{n} u_{i}a_{i} = GCD(a_{1}, …, a_{n}) = d. By theorem 13, we can compute u in time

O(n(log n + log L)M_{B}(L)).  (25)

Note that GCD(u_{1}, …, u_{n}) = 1. As in the proof of lemma 15, we recursively compute the (n − 1) × (n − 1) matrix P′ with det P′ = d′ with first column (u_{1}, …, u_{n−1})^{T}. Then we compute s, t with sd′ + tu_{n} = d, in time M_{B}(L) log L. Finally, we compute the last row of P in time O(nM_{B}(L)). Hence if T (n, L) is the complexity of computing P

T (n, L) = T (n − 1, L) + O((n + log L)M_{B}(L)).

Hence T (n, L) = O(n(n + log L)M_{B}(L)), and this dominates the complexity in (25).

(ii) Given a and P, a straightforward multiplication gives bP in time O(n^{2}M_{B}(L)). But a better bound is needed. For i = 2, …, n, the ith column of P has the form

p_{i} = (\frac{u_{1}}{r_{i}} , \frac{u_{2}}{r_{i}} , …, \frac{u_{i−1}}{r_{i}} , s_{i}, 0, …, 0)^{T}

where r_{i}, s_{i} are the elements described in the proof of lemma 15. So

bP = (⟨b, p_{1}⟩, ⟨b, p_{2}⟩, …, ⟨b, p_{n}⟩)

where ⟨b, p_{i}⟩ indicates scalar product. We will compute the entries of bP in a particular order. Notice that

log |⟨b, p_{i}⟩| = O(L + log n).  (26)

Clearly ⟨b, p_{2}⟩ can be computed in time O(M_{B}(L)). For i = 2, …, n − 1, ⟨b, p_{i+1}⟩ can be obtained from ⟨b, p_{i}⟩ in time O(M_{B}(L + log n)) using the formula

⟨b, p_{i+1}⟩ = \frac{(⟨b, p_{i}⟩ − b_{i}s_{i})r_{i} + b_{i}u_{i}}{r_{i+1}} + b_{i+1}s_{i+1}.

Finally, the first entry ⟨b, p_{1}⟩ in bP can be computed from the last entry in time O(M_{B}(L + log n)). The entire computation costs O(nM_{B}(L + log n)) as claimed. Q.E.D.

Continuing in our pursuit of the matrix U, we now need to compute the matrix P^{∗} such that (aP)P^{∗} will have zero everywhere except the first entry. Clearly, if aP = (g_{1}, g_{2}, g_{3}, …, g_{n}) then g_{1} = d and

P^{∗} = [ \begin{matrix}1 & −g_{2}/d & −g_{3}/d & ⋯ & −g_{n−1}/d & −g_{n}/d \\ 0 & 1 & 0 & ⋯ & 0 & 0 \\ ⋮ & & & ... & & ⋮ \\ 0 & 0 & 0 & ⋯ & 1 & 0 \\ 0 & 0 & 0 & ⋯ & 0 & 1 \end{matrix} ] .

Note that P^{∗} is integer unimodular since d divides each g_{i}. The entries in the first row of P^{∗} has bit-size of O(L + log n). We return to the main problem (20):

Lemma 17 (Hermite Reduction Step) Let a ≠ 0 be the first row of A ∈ ℤ^{m×n} and L = lg ∥A∥_{∞}. Define the matrices P, P^{∗}, U relative to a, as before.

(i) We can compute the matrices P, P^{∗}, U from a in time O(n(n + log L)M_{B}(L + log n)).

(ii) If b is any row of A, we can compute bU from b, P, P^{∗} in time O(nM_{B}(L + log n)).

(iii) We can compute AU from A, P, P^{∗} in time O(mnM_{B}(L + log n)).

(iv) We can compute P, P^{∗}, AU from A in time O(n(n + m + log L)M_{B}(L + log m)).

Proof. (i) We use the previous lemma to compute P in time (23). Similarly, we can compute aP = (g_{1}, …, g_{n}) in time (24). Assuming that we only represent the first row of P^{∗} explicitly, we can also construct P^{∗} in time (24). Next, each entry of U = P P^{∗} can be computed in time O(M_{B}(L + log n)) or

O(n^{2}M_{B}(L + log n))

for the entire matrix U. Thus the three matrices P, P^{∗}, U can be computed in time

O(n(n + log L)M_{B}(L)) + nM_{B}(L + log n) + O(n^{2}M_{B}(L + log n)) = O(n(n + log L)M_{B}(L + log n)).

(ii) To compute bU, we first compute bP in time (24). Then compute (bP)P^{∗} within the same time bound (24). (iii) This amounts to repeating part (ii) m times. (iv) This just adds up parts (i) and (iii). Q.E.D.

Application. Although we could use this lemma repeatedly to compute the Hermite normal form, it seems that the best bounds for bit sizes of intermediate values are exponential in min{m, n}. So this application is only useful for m ≫ n or n ≫ m. We describe another application here. Suppose the columns of A = [a_{1}, …, a_{n}] ∈ ℤ^{m×n} are not linearly independent. Consider the problem of computing B = [b_{1}, …, b_{n−1}] ∈ ℤ^{m×(n−1)} such that Λ(A) = Λ(B). For instance, this is useful in preprocessing a basis before calling the LLL-algorithm, since the LLL-algorithm requires the input columns to be linearly independent.

Note that there exists x ∈ ℤ^{n×1} such that Ax = 0. Finding such an x from A is easily reduced to Gram-Schmidt orthogonalization. In the notations of §IX.1, if A^{∗} = [a_{1}^{∗}, …, a_{n}^{∗}] is the Gram-Schmidt version of A, then a_{i}^{∗} = 0 for some i. This means

a_{i}^{∗} = 0 = a_{i} − ∑_{j=1}^{i−1} μ_{ij}a_{j}^{∗}.

Recall that μ_{ij} are rational numbers. It is then easy to find integers t_{1}, …, t_{i} such that

0 = t_{i}a_{i} + ∑_{j=1}^{i−1} t_{j}a_{j}.  (27)

We may set x = (t_{1}, …, t_{i}, 0, …, 0)^{T}. Clearly, we may assume that GCD(x) = 1. To find B, we can use the Hermite reduction step to give us a unimodular matrix U ∈ ℤ^{n×n} such that Ux = (1, 0, …, 0)^{T}. Since (AU^{−1})(Ux) = 0, we conclude that the first column of AU^{−1} is zero. The desired matrix B may be taken to comprise all but the first column of AU^{−1}.

Algebraic complexity of the Hermite Reduction Step.

We have given an upper bound on the bit complexity of the Hermite Reduction Step (20). But suppose we want its complexity in an algebraic model of computation (§0.6). It is clear from the preceding that problem (20) over ℤ can be solved in polynomial time in an algebraic model M if the elementary operations or basis of M includes the following:

+, −, , ×, cGCD(x, y).  (28)

Here cGCD denotes the co-GCD primitive (§II.2) that returns a co-factor (s, t) of an input pair (x, y) ∈ ℤ^{2}: sx + ty = GCD(x, y). Hafner and McCurley [8, p. 1075] suggested that there may be no solution in case our basis comprises only the ring operations (+, −, ×), i.e., if we omit cGCD. Let us show this. If (20) can be solved in the algebraic model M then GCD(x, y) can be solved in M in constant time by reduction to the n = 2 case. This assumes (as we may) that the ring operations are available in M. Suppose π(x, y) is a branching program in M for computing the GCD(x, y). The inputs x, y and constants used in π are all integers, and at each decision node, there is an integer polynomial f(X, Y) whose sign at the input (x, y) determines the flow of computation. The finiteness of π means that there are finitely many leaf nodes in the branching program. At each leaf ℓ, there is an integer polynomial P_{ℓ}(X, Y) such that if input (x, y) terminates at ℓ the P_{ℓ}(x, y) = GCD(x, y). Let S_{ℓ} denote the set of all (x, y) ∈ ℝ^{2} that terminates at ℓ. Note that it makes sense to feed pairs (x, y) of real numbers to π and hence the set S_{ℓ} is well-defined. Clearly, S_{ℓ} is a semi-algebraic set (i.e., defined by a finite set of polynomial inequalities). By basic properties of semi-algebraic sets, there is some leaf ℓ such that S_{ℓ} has the following properties: S_{ℓ} contains an infinite cone C which in turn contains infinitely many vertical rays of the form R_{i} = {(x_{i}, y) : y ≥ c_{i}} where x_{i} ∈ ℤ is prime and c_{i} is arbitrary, for i = 0, 1, 2, . . .. Focus on the output polynomial P_{ℓ}(X, Y) at such an ℓ. We may pick a ray R_{i} such that such that none of the non-zero coefficients of the Y ’s in P_{ℓ}(X, Y) vanish when X is replaced by x_{i}. Since there are infinitely many prime y’s such that (x_{i}, y) ∈ R_{i} ⊆ S_{ℓ}, we conclude that P_{ℓ}(x_{i}, Y) is the constant 1 and P_{ℓ}(X, Y) does not depend on Y. Next suppose P_{ℓ} has X-degree < d. Pick a disc D ⊆ S_{ℓ} large enough so that there are prime numbers y_{0} and u_{i} (i = 1, …, d) such that D contains (u_{i}, y_{0}) for i = 1, …, d and none of the non-zero coefficients of X in P_{ℓ}(X, Y) vanishes when y_{0} replaces Y. Again, we argue that P_{ℓ}(X, y_{0}) is the constant 1 and P_{ℓ}(X, Y) does not depend on X. Thus P_{ℓ}(X, Y) must be the constant 1. But, clearly S_{ℓ} contains a pair (a, b) of integers such that GCD(a, b) > 1. This is a contradiction since 1 = P_{ℓ}(a, b) = GCD(a, b) > 1.

§7. Bachem-Kannan Algorithm  Lecture X  Page 23

Exercises Exercise 6.1: (i) Show that the matrix P in lemma 15 is a product of elementary integer unimod-

ular matrices (§VIII.1). HINT: use induction on n. (ii) Do the same for P^{∗}. Hence conclude that U in (20) is a product of elementary integer unimodular matrices. □

Exercise 6.2: (i) What is the algebraic complexity of the Hermite reduction step assuming the basis (28)? (ii) Show that if 2×2 Smith normal form (see §8) can be solved in constant time in an algebraic model M then the general Smith normal form problem can be solved in polynomial time in M. (iii) Assume an algebraic computation model M whose basis is given by (28). Show that the 2 × 2 Smith normal form problem is equivalent to the following: given a, b, c ∈ ℤ find U, V, W, Z ∈ ℤ such that aUV + bV W + cW Z = GCD(a, b, c). HINT: first reduce the 2 × 2 matrix to Hermite normal form. (iv) (Open) Prove that the problem in (iii) cannot be solved in constant time in model M. □

Exercise 6.3: In the above application, work out efficient algorithms for:

(i) Computing the integers t_{1}, …, t_{n} in (27) from μ_{ij}’s. (ii) Computing the inverse U^{−1} from U. □

Exercise 6.4: Suppose a = (a_{1}, …, a_{n})^{T} and b = (b_{1}, …, b_{n})^{T} (n ≥ 2) be two columns. Under what circumstances is there a unimodular matrix U whose first and second columns are a and b? A necessary condition is that GCD(a) = GCD(b) = 1 and the GCD of all minors of order 2 of [a, b] is 1. □

Exercise 6.5: (Bass) Let M = [a_{ij}] be an n × n matrix over a ring R. Say M is invertible iff its determinant is invertible in R. But det M = a_{11}A_{11} + ... + a_{1n}A_{1n} where A_{ij} is the co-factor of a_{ij}. Write a_{i} for a_{1i}. Then invertibility of M implies R = Ideal(a_{1}, …, a_{n}). We also call (a_{1}, …, a_{n}) a unimodular row in this case. The converse asks: is every unimodular row the row of an invertible matrix? This is related to a conjecture of Serre’s which was answered affirmatively by Quillen and Suslin. Verify the counterexample: R = ℝ[X, Y, Z]/(X^{2} + Y^{2} + Z^{2} = 1). The unimodular row (\overline{X}, \overline{Y} , \overline{Z}) is not the first row of any invertible matrix. Here \overline{u} denote the image of u in the canonical map from ℝ[X, Y, Z] to R. □

§7. Bachem-Kannan Algorithm

We present a polynomial time algorithm for Hermite normal form. It is essentially that of Bachem and Kannan [1], extended to non-square matrices of arbitrary rank. We assume that the input matrix A ∈ ℤ^{m×n} has been preprocessed so that the ith principal minor is non-zero for i = 1, …, ρ where ρ > 0 is the rank of A. The Extended Bareiss Algorithm (§2) can convert any matrix A into this form. Moreover, the complexity of this conversion is dominated by the subsequent computation.

Algorithm. The algorithm is relatively straightforward to describe. For any matrix M, let

⟨M⟩_{i}:= M(1, 2, …, i; 1, 2, …, i)

using the general notation of §I.10. So ⟨M⟩_{i} is just the ith principal submatrix of M. On input A, th algorithm has n stages where in the sth (s = 1, …, n) stage, we seek to put ⟨A⟩_{s} into the Hermite normal form. Stage 1 requires no action. Suppose that we have just completed the (s−1)st stage. The sth stage has two phases. Elimination phase: we eliminate (i.e., zero out) the first s − 1 entries in the sth column. This will make ⟨A⟩_{s} a lower triangular matrix. Reduction phase: we then reduce the off-diagonal entries of ⟨A⟩_{s} so that each such entry is non-negative and less than the corresponding diagonal entry in its row. This completes the sth stage.

Bachem-Kannan Algorithm

Input: A ∈ ℤ^{m×n} and ρ ≥ 1 the rank of A.

Assume the ith principal minor is non-zero for i = 1, …, ρ

Output: H ∈ ℤ^{m×n}, the HNF of A, and U ∈ ℤ^{n×n}, a unimodular matrix

such that H = AU.

INITIALIZATION:

1. H:= A; U := I the identity matrix. MAIN LOOP:

2. for s ← 2 to n do ELIMINATION PHASE:

3. for i ← 1 to min{s − 1, ρ} do

4. By postmultiplying H with a suitable unimodular matrix K,

eliminate the (i, s)th entry of H; Update U ← UK.

REDUCTION PHASE:

5. if s > ρ, skip this phase; else continue.

6. for j ← s − 1 downto 1 do

7. for i ← j + 1 to s do

8. By postmultiplying H with a suitable unimodular matrix K, reduce the (i, j)th entry of H; Update U ← UK.

Call step 4 an “elimination step” and step 8 a “reduction step”. Note that when s > ρ, the reduction phase is omitted since column s would be zero after the elimination phase. The order of reduction represented by the double for-loop (steps 6 and 7) is important for the analysis: this is an improvement from Chou and Collins [4].

The reduction step is rather obvious: the entry x to be reduced is replaced by x mod d where d is the diagonal element in the same row. The rest of the column of x is modified accordingly. We now illustrate the elimination step: it is basically the 2 × 2 version of the Hermite reduction step (§6). For instance, suppose H ∈ ℤ^{5×6} already has its 3rd principal submatrix in HNF. The following shows the first elimination step of stage 4:

H = [h_{1}, h_{2}, h_{3}, h_{4}, h_{5}, h_{6}] → [h′_{1}, h_{2}, h_{3}, h′_{4}, h_{5}, h_{6}] = H′,

[ \begin{matrix}a_{1,1} & 0 & 0 & a_{1,4} & ∗ & ∗ \\ a_{2,1} & a_{2,2} & 0 & a_{2,4} & ∗ & ∗ \\ a_{3,1} & a_{3,2} & a_{3,3} & a_{3,4} & ∗ & ∗ \\ a_{4,1} & a_{4,2} & a_{4,3} & a_{4,4} & ∗ & ∗ \\ ∗ & ∗ & ∗ & ∗ & ∗ & ∗ \end{matrix} ] → [ \begin{matrix}a′_{1,1} & 0 & 0 & 0 & ∗ & ∗ \\ a′_{2,1} & a_{2,2} & 0 & a′_{2,4} & ∗ & ∗ \\ a′_{3,1} & a_{3,2} & a_{3,3} & a′_{2,4} & ∗ & ∗ \\ a′_{4,1} & a_{4,2} & a_{4,3} & a′_{2,4} & ∗ & ∗ \\ ∗ & ∗ & ∗ & ∗ & ∗ & ∗ \end{matrix} ] .

As in (21) (§6), this amounts to replacing columns h_{1}, h_{4} (respectively) by h′_{1}, h′_{4} which is defined as follows:

[h′_{1}, h′_{4}] ← [h_{1}, h_{4}] [ \begin{matrix}s & a_{1,4}/g \\ t & −a_{1,1}/g \end{matrix} ]

where g = GCD(a_{1,1}, a_{1,4}) = s · a_{1,1} + t · a_{1,4}. We may assume that (see §5)

|s| < |a_{1,4}|,  |t| ≤ |a_{1,1}|.  (29)

We may write this elimination step as H′ = HK where

K = [ \begin{matrix}s & 0 & 0 & a_{1,4}/g & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ t & 0 & 0 & −a_{1,1}/g & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{matrix} ] .  (30)

If U is the unimodular matrix such that H = AU then we update U via U ← UK.

Analysis. There are subtleties in proving a polynomial running time. In this analysis, we simply write ∥M∥ instead of ∥M∥_{∞} for the ∞-norm of a matrix M. With

λ_{0} := ∥A∥, L:= lg λ_{0}.

The analysis amounts to showing that throughout the algorithm, lg ∥H∥ and lg ∥U∥ are polynomiallybounded in terms of m, n and L.

Bound between stages. Let H^{(s)} denote the H-matrix just after the sth stage. Thus H^{(1)} is equal to the initial matrix A. Also let the unimodular matrix that transforms A to H^{(s)} be denoted U^{(s)}:

H^{(s)} = AU^{(s)}.

So U^{(1)} is the identity matrix. Letting

λ_{1} :=(ρλ_{0})^{ρ},  (31)

we see that every minor of A is bounded by λ_{1}. From this we infer

∥⟨H^{(s)}⟩_{s}∥ ≤ λ_{1}.  (32)

This is because each step (elimination or reduction) in the first s stages does not change the sth principal minor of H, and H is initially equal to A. Since ⟨H^{(s)}⟩_{s} is lower triangular, this means that each diagonal entry of ⟨H^{(s)}⟩_{s} is bounded by λ_{1}. But the off-diagonal entries are also bounded by λ_{1}, since ⟨H^{(s)}⟩_{s} is in HNF. Thus (32) is verified.

First assume s ≤ ρ. Note that U^{(s)} has the form

U^{(s)} = [ \begin{matrix}⟨U^{(s)}⟩_{s} & 0 \\ 0 & I \end{matrix} ] .  (33)

So

⟨U^{(s)}⟩_{s} = ⟨A⟩_{s}^{−1}⟨H^{(s)}⟩_{s}.  (34)

Each entry of ⟨A⟩_{s}^{−1} is bounded by λ_{1}/ det(⟨A⟩_{s}) (see §1). So a typical element U_{ij}^{(s)} of ⟨U^{(s)}⟩_{s} is bounded by

|U_{ij}^{(s)}| ≤ ∑_{k=1}^{s} \frac{λ_{1}}{| det⟨A⟩_{s}|}|H_{kj}^{(s)}| ≤ λ_{1},

since

∑_{k=1}^{s} |H_{kj}^{(s)}| ≤ 1 + ∑_{k=j}^{s}(|H_{kk}^{(s)}| − 1)  (35)

≤ ∏_{k=1}^{s} |H_{kk}^{(s)}|  (36)
= | det⟨A⟩_{s}|.

Inequality (35) is a consequence of Hermite normal form of ⟨H^{(s)}⟩_{s}. Inequality (36) exploits the elementary bound

1 + ∑_{k=1}^{ℓ}(x_{k} − 1) ≤ ∏_{k=1}^{ℓ} x_{k}  (37)

which holds for any ℓ ≥ 1 positive integers x_{1}, …, x_{ℓ} (see Exercise). From (33), we infer that

∥U^{(s)}∥ ≤ ∥⟨U^{(s)}⟩_{s}∥ ≤ λ_{1}.  (38)

From H^{(s)} = AU^{(s)} and the special form (33), we conclude

∥H^{(s)}∥ ≤ λ_{2}  (39)

where

λ_{2} := ρλ_{0}λ_{1} = (ρλ_{0})^{1+ρ}.  (40)

We have therefore given a bound on ∥U^{(s)}∥ and ∥H^{(s)}∥ for s ≤ ρ.

Now let s > ρ. Clearly H^{(s)} still satisfies the bound (39) since in stage s we eliminated column s while the remaining columns are unchanged from H^{(s−1)}. We conclude that ∥H∥ ≤ λ_{2} holds in the transition between any two successive stages.

What about U^{(s)}? Eliminating the (r, s)th entry amounts to adding some multiple c_{r,s} of column r (r ≤ ρ) to column s. The multiple c_{r,s} is bounded by λ_{0}. Therefore, it increases ∥U^{(s)}∥ by a factor of (λ_{0} + 1). We perform ρ such elimination steps to entirely eliminate column s. Thus

∥U^{(s)}∥ ≤ (λ_{0} + 1)^{ρ}∥U^{(s−1)}∥.

The maximum size bound is when s = n:

∥U^{(n)}∥ ≤ (λ_{0} + 1)^{ρ(n−ρ)}∥U^{(ρ)}∥ ≤ λ_{3},  (41)

where

λ_{3} :=(λ_{0} + 1)^{ρ(n−ρ)}λ_{1}.  (42)

Bounds on H during a stage. What remains is to bound ∥H∥ and ∥U∥ during a stage. In other words, although the entries in H and U are nicely bounded between two successive stages, we must ensure that they do not swell up excessively within a stage. In our analysis, we shall use the observation that once the (s, s)th element is “fixed” at the end of the sth stage, it is hereafter bounded by λ_{1}. If it changes at all, it can only become smaller (replaced by a divisor). In fact, the product of all the “fixed” diagonal elements is bounded by λ_{1}. Let us now focus on stage s for some s = 1, …, n. Of course, the columns of interest are the first min{s − 1, ρ} columns and column s.

ELIMINATION PHASE: let H^{(r,s)} be the H matrix just before the (r, s)th entry is eliminated (r = 1, …, min{s − 1, ρ}). Let h_{j}^{(r,s)} be the jth column of H^{(r,s)} and H_{i,j}^{(r,s)} be the (i, j)th entry of H^{(r,s)}. Initially, we have

∥h_{j}^{(1,s)}∥ ≤ \{ \begin{matrix}λ_{2}, & for j = 1, …, min\{s − 1, ρ\}, \\ λ_{0}, & for j = s. \end{matrix}

When we eliminate the (r, s)th entry, column s is transformed (cf. equation (30)) according to the rule

h_{s}^{(r+1,s)} ← \frac{H_{r,s}^{(r,s)}h_{r}^{(r,s)} − H_{r,r}^{(r,s)}h_{s}^{(r,s)}}{g}  (43)

where g = GCD(H_{r,r}^{(r,s)}, H_{r,s}^{(r,s)}). At the same time, column r is transformed according to the rule

h_{r}^{(r+1,s)} ← c · h_{r}^{(r,s)} − c′ · h_{s}^{(r,s)}  (44)

where |c| < |H_{r,s}^{(r,s)}| and |c′| ≤ |H_{r,r}^{(r,s)}| (cf. equation (29)). Define

β_{r}:= λ_{2}(λ_{1} + λ_{2})^{r−1}.

Inductively assume that column s is bounded as follows:

∥h_{s}^{(r,s)}∥ ≤ β_{r}.  (45)

This is true for r = 1. From (43), we extend the inductive hypothesis to r + 1:

∥h_{s}^{(r+1,s)}∥ ≤ |H_{r,s}^{(r,s)}| · ∥h_{r}^{(r,s)}∥ + |H_{r,r}^{(r,s)}| · ∥h_{s}^{(r,s)}∥
≤ β_{r} · λ_{2} + λ_{1} · β_{r}
= β_{r}(λ_{1} + λ_{2}) = β_{r+1}.

From (44), we similarly obtain a bound on column s:

∥h_{r}^{(r,s)}∥ ≤ β_{r}.  (46)

Let H^{(s,s)} be the H matrix just after the (s − 1, s)th entry of H is eliminated. Then the bounds (45) and (46) extend to

∥h_{s−1}^{(s,s)}∥ ≤ β_{s}, ∥h_{s}^{(s,s)}∥ ≤ β_{s}.

We conclude that throughout an elimination phase, each entry of H is bounded by

β_{s} = λ_{2}(λ_{1} + λ_{2})^{s−1} < λ_{4}

where

λ_{4} :=(2λ_{2})^{ρ} = 2^{ρ}(ρλ_{0})^{ρ(1+ρ)}.  (47)

REDUCTION PHASE: So H^{(s,s)} is the H matrix just before the start of the reduction phase. Note that we may assume s ≤ ρ. Let h_{j}^{(s)} be the jth column of H^{(s,s)}. Also let ĥ_{j}^{(s)} be the jth column at the end of the reduction phase: these are called the “reduced vectors”. Note that the reduced vectors are columns of H^{(s)} and hence satisfy

∥ĥ_{j}^{(s)}∥ ≤ λ_{2}.  (48)

Exploiting the special sequencing of reduction steps in the algorithm (following Chou-Collins), we see that

ĥ_{j}^{(s)} = h_{j}^{(s)} − ∑_{r=j+1}^{s} b_{r,j,s}ĥ_{r}^{(s)}  (49)

for suitable constants b_{r,j,s}. The point is that reduced vectors appear on the right-hand side of (49). The entries of column j in H is bounded by λ_{4} at the start of the reduction phase. To reduce column j (j = 1, …, s−1), we first reduce its j +1st entry is reduced by subtracting the column b_{j+1,j,s}ĥ_{j+1}^{(s)} (cf. equation (49)). Clearly b_{j+1,j,s} ≤ λ_{4} so that entries of column j are bounded by λ_{4}(1 + λ_{2}) after this reduction step. Inductively, it is easy to see that after the j + kth (k = 1, 2, …) entry of column j is reduced, the entries of column j are bounded by

λ_{4}(1 + λ_{2})^{k}.

So just after the s − 1st entry is reduced, its entries are bounded by

λ_{4}(1 + λ_{2})^{s−1−j} < λ_{4}(2λ_{2})^{ρ} = λ_{4}^{2}.

Finally, we reduce the sth entry of column j. But the result is a reduced column bounded as in (48). It follows that the bound

∥H∥ ≤ λ_{4}^{2}

holds throughout the stage.

§8. Smith Normal Form  Lecture X  Page 28

Bounds on U during a stage. First assume s ≤ ρ. It suffices to use the bound

∥U∥ ≤ n∥A^{−1}∥ · ∥H∥ ≤ nλ_{1}λ_{4}^{2}

since the relation U = A^{−1}H holds throughout the stage. holds throughout the stage. Suppose s > ρ. There is no reduction phase and the argument showing ∥U^{(s)}∥ ≤ λ_{3} in (41) actually shows that ∥U∥ ≤ λ_{3} throughout stage s. This concludes our analysis.

We summarize the foregoing analysis: the entries of H and U matrices are uniformly bounded by

λ_{5} := λ_{3} + λ_{4} < λ_{1}((2λ_{0})^{ρ(n−ρ)} + nλ_{4}^{2})

throughout the algorithm. Since L = lg λ_{0}, we get

L′:= lg λ_{5} = O(ρnL + ρ^{2} lg ρ)

as a bound on the bit-size of entries. There are O(ρ^{2}n) column operations and O(ρ^{2}) co-GCD operations on matrix entries. Each column operation takes O(m) ring operations on matrix entries. Hence the cost of co-GCD operations is dominated by the cost of column operations, which amounts to:

Theorem 18 With L = lg λ_{0}, the matrix entries have bit-size that are uniformly bounded by

L′ = O(ρ[nL + ρ lg ρ)])

in the Bachem-Kannan algorithm. The bit-complexity of the algorithm is O(mnρ^{2}M_{B}(L′)).

Remarks. Our complexity analysis is somewhat more involved than Bachem and Kannan [1], in part because the input matrix is non-square and has dependent rows and columns. For instance, if n = ρ then lg(λ_{3}) = O(n(L + lg n) and not O(n^{2}L).

Suppose we want to compute the HNF H(A) of an arbitrary matrix A. Accordingly, we submit A to the Extended Bareiss algorithm (§2) to obtain B = P AQ where P and Q are permutation matrices. Now we submit B to the Bachem-Kannan algorithm which output H = H(B). computes the HNF of B. We leave it as an exercise to show:

H(A) = P^{−1}H.  (50)

Exercises

Exercise 7.1: (Chou-Collins) Verify the inequality (37) □

Exercise 7.2: Show the equation (50). HINT: this depends on the particular structure of our Extended Bareiss algorithm. □

Exercise 7.3: Analyze the algorithm assuming A is square and non-singular. □

§8. Smith Normal Form

H.J.S. Smith (1861) introduced the normal form bearing his name. Let A ∈ ℤ^{m×n}. We say A is in Smith normal form (SNF) if it is diagonal, that is, (A)_{i,j} = 0 for i ≠ j, and its diagonal entries are non-negative satisfying

(A)_{i−1,i−1}|(A)_{i,i},  for i = 2, …, min{m, n}.  (51)

Since every number divides 0, but 0 divides only itself, we conclude from (51) that if (A)_{i,i} = 0 for some i, then (A)_{j,j} = 0 for all j > i. The multi-set of non-zero diagonal entries of a Smith normal form matrix A,

{(A)_{1,1}, (A)_{2,2}, …, (A)_{r,r}}

where (A)_{r,r} > 0 and (A)_{r+1,r+1} = 0, is called the set of Smith invariants of A. We also call (A)_{i,i} the ith Smith invariant (i = 1, …, r). [In the literature, the Smith invariants of A are also called “invariant factors” of A.]

By elementary operations in this section, we mean elementary integer row or column operations. We say that two matrices are equivalent if they are inter-transformable by elementary operations.

Lemma 19 Every integer matrix A can be brought into a Smith normal form S by a sequence of elementary operations.

We leave the proof to an Exercise. The algorithm below implies this, of course, but it is instructive for the student to give a direct proof. We will show that S is unique for A, and so S is the Smith normal form of A, usually denoted S(A). For k = 1, …, min{m, n}, let δ_{k}(A) denote the GCD of all the order k minors of A. In particular, δ_{1}(A) is the GCD of every entry of A.

Lemma 20

(i) The elementary operations on a matrix A preserves δ_{k}(A).

(ii) The set of Smith invariants of A is unique.

(iii) The Smith normal form of A is unique.

Proof. (i) This is immediate from our treatment the γ-invariants for HNF. (ii) Let the rank of A be r. Then δ_{k}(A) ≠ 0 iff k = 1, …, r. From the definition of δ’s, it is clear that

δ_{k}(A)|δ_{k+1}(A)

for k = 1, …, r − 1. Let S be a Smith normal form of A. Since A and S have the same rank, we conclude that (S)_{k,k} ≠ 0 iff k = 1, …, r. Note that δ_{k}(S) = (S)_{1,1}(S)_{2,2} ⋯ (S)_{k,k}. In view of part (i), we conclude

(S)_{1,1}(S)_{2,2} ⋯ (S)_{k,k} = δ_{k}(A).

It follows that (S)_{1,1} = δ_{1}(A) and (S)_{k+1,k+1} = δ_{k+1}(A)/δ_{k}(A) (k = 2, …, r − 1). (iii) This follows from (ii) since the Smith normal form is determined by the set of its Smith invariants.

Q.E.D.

Polynomial-time Algorithm. Computing the Smith normal form of a matrix A is equivalent to computing the set of Smith invariants. For some applications, it is desirable to also know the unimodular matrices U, V such that

S(A) = UAV.

For instance, it is easy to compute the 1st Smith invariant – it is just the GCD of all the matrix entries. But computing U and V such that (UAV)_{1,1} has this invariant is non-trivial. As usual, the difficulty in giving a polynomial time algorithm is the possibility of exponential size intermediate entries. We describe an algorithm from Bachem-Kannan [1], based on a Hermite normal form algorithm.

It suffices to show how to perform a Smith Reduction Step (in analogy to the Hermite Reduction Step of §6): given a matrix A ∈ ℤ^{m×n}, compute two unimodular matrices U ∈ ℤ^{m×m}, V ∈ ℤ^{n×n} such that UAV is “Smith-reduced”. In general, a matrix M is said to be Smith-reduced if:

(i) The first row and first column of M are zero except for the top-left corner entry (M)_{1,1}. (ii) (M)_{1,1} divides all the remaining entries of M.

Our Hermite normal form is based on column operations. We now need the “row version” of the normal form: a matrix A is in Row Hermite normal form (abbr. RHNF) if its transpose A^{T} is in Hermite normal form. Using elementary row operations, or multiplication by unimodular matrix on the left, we can transform any matrix into its RHNF. Algorithms for RHNF is trivially obtained from HNF algorithms, simply by interchanging the roles of rows and columns.

Smith Reduction Step Input: A ∈ ℤ^{m×n}. Assume (A)_{1,1} ≠ 0. Output: Unimodular matrices U ∈ ℤ^{m×m} and V ∈ ℤ^{n×n}, and matrix S

such that S is Smith-reduced and HNF, and S = UAV.

INITIALIZATION: 1. S ← A.

2. U ← I_{m}; V ← I_{n} (identity matrices). MAIN LOOP:

{Loop Invariant: S = UAV } 3. while S is not Smith-reduced do

4. if (S)_{1,1} is the only non-zero entry in the first row and first column, then (S)_{1,1} does not divide some (S)_{i,j}.

5. In this case, add column j to column 1, and update V.

6. Apply a RHNF algorithm to S and update U accordingly.

7. Apply a HNF algorithm to S and update V accordingly.

Analysis. The bit-size of entries remain bounded between successive while-iterations: this is because S is in HNF on exit from an iteration, and so the largest entry lies on the diagonal. But the diagonal entries of S are bounded by the determinant of the input matrix A, since the elementary operations preserve δ_{k}(A) for each k. Since the RHNF and HNF algorithms are polynomial time, we conclude that each while-iteration is polynomial-time.

It remains to show that the number of iterations is polynomial. Note that whenever (S)_{1,1} changes after an iteration, it is being replaced by a factor of itself. We claim that (S)_{1,1} must change at least in every other iteration. To see this, consider the two possibilities: either step 5 (adding column j to column 1) is executed or it is not. Step 5 is executed: then column 1 contains an entry not divisible by (S)_{1,1} before the RHNF algorithm. After the RHNF algorithm, (S)_{1,1} will contain the the GCD of entries in column 1, and this will be a proper factor of its previous value. Step 5 is not executed: then row 1 or column 1 must have at least some other non-zero entry. Again there are two cases. (i) If all of these non-zero entries are divisible by (S)_{1,1} then after this iteration, (S)_{1,1} is the only non-zero entry in row 1 and column 1. If this is not the last iteration, then we already saw that (S)_{1,1} will change in the next iteration. (ii) If there is an entry that is not divisible by (S)_{1,1} then either the RHNF or HNF algorithm will change (S)_{1,1} to a smaller value. Hence the number of iterations is at most 1 + 2 lg ∥A∥.

One other remark: the output matrices U, V is polynomial size. This is because each matrix is the product of O(lg ∥A∥) many component matrices (produced by the call to HNF or RHNF or by Step 5). But each component matrix is polynomial size. This concludes our analysis of the Smith Reduction Step algorithm. Clearly the Smith normal form of an m × n matrix can be reduced to min{m, n} Smith Reduction steps. Moreover, the result of Smith Reduction Step is HNF, and hence polynomially bounded in terms of the original input. This proves:

Theorem 21 There is a polynomial-time algorithm to compute the Smith normal form S(A) of a matrix A and also unimodular matrices such that S(A) = UAV.

Exercises

Exercise 8.1: Show lemma 19 by a direct argument (i.e., without reference to the existence of the SNF algorithm). □

Exercise 8.2: Let d_{i} (i = 1, …, r) be the ith Smith invariant of A. Write

d_{i} = p_{1}^{e_{i,1}}p_{2}^{e_{i,2}}, …, p_{ℓ_{i}}^{e_{i,ℓ_{i}}} , (ℓ_{i} ≥ 1)

where p_{i} is the ith prime number. Call the prime power p_{j}^{e_{i,j}} an elementary divisor of A. Show that two matrices are equivalent iff they have the same rank and the same set of elementary divisors. □

Exercise 8.3: Analyze the complexity of the Smith Reduction Step and the associated SNF algorithm. □

Exercise 8.4: Let A, B, C ∈ ℤ^{n×n}. We say A is irreducible if, whenever A = BC then either B or C is unimodular. Otherwise A is reducible. If A = BC, we call C a right divisor of A or C right-divides A. Write B|A in this case. Similarly, there is a notion of left divisor.

(i) An irreducible matrix is equivalent to diag(1, …, 1, p), the diagonal matrix whose main diagonal has all ones except the last entry is a prime p. (ii) A necessary and sufficient condition for a square matrix to be irreducible is that its determinant is prime. (iii) A reducible matrix can be factored into a finite product of irreducible matrices. Formulate a uniqueness property for this factorization. □

Exercise 8.5: Let A, B, C, D ∈ ℤ^{n×n}. Assume A and B are not both zero. Call D a (right) greatest common divisor (GCD) of A, B if D|A and D|B and for any other C that divides both A and B, then C|D. (See definitions in previous exercise.)

(i) Show that GCDs of A, B exist. Moreover, if D is a GCD then D = P A + QB for some P, Q. HINT: consider the SNF of [A|B]. (ii) If D, D′ are two GCD’s of A then D = UD′ for some unimodular matrix U. □

§9. Further Applications  Lecture X  Page 32

§9. Further Applications

Linear Diophantine equations. Let A ∈ ℤ^{m×n} and b = (b_{1}, …, b_{n}) ∈ ℤ^{1×n}. Consider the problem of solving the linear system

x · A = b  (52)

for an unknown x = (x_{1}, …, x_{m}) ∈ ℤ^{1×m}. Such a system is also called a Diophantine linear system. For simplicity, assume m ≥ n; otherwise, the n equations in (52) are not independent and some may be omitted. Let S = UAV be the Smith normal form of A and let the diagonal elements of S be d_{1}, …, d_{n}. Then (52) implies

(xU^{−1})(UAV) = bV,
x̂S = b̂.

where x̂ = (x̂_{1}, …, x̂_{m}) = xU^{−1} and b̂ = (b̂_{1}, …, b̂_{n}) = bV. The last equation amounts to

x̂_{i}d_{i} = b̂_{i},  i = 1, …, n.  (53)

Suppose d_{1}, …, d_{r} are non-zero and d_{r+1} = ⋯ = d_{n} = 0. Then the system (52) has solution iff (i) d_{i}|b̂_{i} for i = 1, …, r and (ii) b̂_{i} = 0 for i = r + 1, …, n. If these conditions are satisfied then a general solution to (53) is given by setting x̂_{i} = b̂_{i}/d_{i} for i = 1, …, r and arbitrary assignments to x̂_{i} for i = r + 1, …, m. For instance, we may choose

x̂ = (b̂_{1}/d_{1}, …, b̂_{r}/d_{r}, 0, …, 0).

From any such particular solution, we obtain a solution x = x̂U to the original system (52).

Homogeneous Case. The important special case of (52) where b = 0 is said to be homogeneous. A solution to this homogeneous system xA = 0 is called a null-vector of A. The set N(A) ⊆ ℤ^{m} of these null-vectors forms a ℤ-module. By the Hilbert basis theorem for modules (see Lecture XI), N(A) has a finite basis. Now e_{r+1}, …, e_{m} is a basis for the set of solutions to (53) in the homogeneous case. Here e_{i} is the m-vector with 1 in the ith position and 0 everywhere else. We conclude that the set

e_{r+1}U, e_{r+2}U, …, e_{m}U

is a basis for N(A). Therefore, we may consider a complete solution of (52) to have the form (x̂U, e_{r+1}U, …, e_{m}U) where x̂ is an particular solution to (53). Notice that N(A) is also a lattice, albeit represented in “dual form”.

If A ∈ ℤ^{m×n} and m > n then we have just shown that the homogeneous diophantine system x · A = 0 has non-trivial solutions. An interesting question is whether there exists small integer solutions. Siegel (1929) shows: there is an x ∈ N(A) satisfying

∥x∥ < 1 + (m∥A∥)^{n/(m−n)}.  (54)

We leave the demonstration to an exercise.

Finitely-generated Abelian groups. Smith normal forms are intimately related to the theory of finitely generated Abelian groups. Let G be a finitely-presented Abelian group, that is, G is represented by n generators x_{1}, …, x_{n} and m relations of the form

∑_{j=1}^{n} a_{ij}x_{j} = 0, (a_{ij} ∈ ℤ)

for i = 1, …, m. (Note the convention of writing the operations of an Abelian group “additively”.) The corresponding relations matrix is an m × n integer matrix A where a_{ij} = (A)_{i,j}. We may rewrite the relations in the form Ax = 0 where x = (x_{1}, …, x_{n})^{T}. In the special case where m = 0 (alternatively, the matrix A is all zero) deserves mention: the group G in this case is called the free Abelian group of rank n. Clearly, G ≈ ℤ^{n} where ≈ indicates group isomorphism.

Let the Smith normal form of A be S = S(A) = UAV for some unimodular matrices U, V. This amounts to transforming the generators of G to (y_{1}, …, y_{n})^{T} = V^{−1}(x_{1}, …, x_{n})^{T}. Then S · (y_{1}, …, y_{n})^{T} = 0. If S has rank r and the diagonal elements of S are d_{1}, …, d_{min(m,n)} then we see that d_{i}y_{i} = 0 for i = 1, …, r and y_{r+1}, …, y_{n} satisfies no relations whatsoever. Each y_{i} generates the subgroup G_{i} = ℤy_{i} of G. Moreover, G_{i} ≈ ℤ_{d_{i}} for i = 1, …, r and G_{i} ≈ ℤ for i = r + 1, …, n. Clearly G is a direct sum of the G_{i}’s: G = ⊕_{i=1}^{n}G_{i}. There are three kinds of subgroups: d_{i} = 0: these corresponds to torsion-free subgroups G_{i} ≈ ℤ. The number β of these subgroups is called the Betti number of G. d_{i} = 1: these are trivial subgroups, and may be omitted in the direct sum expression. d_{i} ≥ 2: these gives rise to finite cyclic groups G_{i}. These d_{i}’s are called torsion coefficients of G.

We have just proven the “fundamental theorem of finitely generated Abelian groups”: every finitely presented Abelian group G on n generators can be written as a direct sum G = ⊕_{i=0}^{r}H_{i} where H_{0} is a free Abelian group of rank β, and each H_{i} (i = 1, …, r) is a finite cyclic group of order d_{i} ≥ 2 satisfying such that d_{1}|d_{2}| ⋯ |d_{r}. The numbers β, d_{1}, …, d_{r} are uniquely determined by G.

Exercises

Exercise 9.1: Let a, b, c ∈ ℤ where a, b are relatively prime. Suppose sa + tb = 1. Show that the general solution of the Diophantine equation ax + by = c is (x, y) = (sc + nb, tc − na) where n is any parameter. □

Exercise 9.2: Consider N(A) in case n = 1. Say A = (a_{1}, …, a_{m})^{T}.

(i) Let s = (s_{1}, …, s_{m}) be a co-factor of (a_{1}, …, a_{m}). Show that ℤ · N(A) + ℤ · s is the unit lattice ℤ^{m}. (ii) For 1 ≤ i < j ≤ m let T (i, j) be the m-vector that is zero everywhere except a_{i}/GCD(a_{i}, a_{j}) at the jth position and −a_{j}/gcd(a_{i}, a_{j}) at the ith position. The set of these T (i, j)’s generates N(A). □

Exercise 9.3: (i) Show the bound of Siegel (54). HINT: for H a parameter to be chosen, let C = C_{H} be the cube comprising points x ∈ ℝ^{m} where ∥x∥ ≤ H. Let α : ℝ^{m} → ℝ^{n} be the linear map given by α(x) = x · A. Give a cube C′ in ℝ^{n} the contains α(C). Using a pigeon hold argument to show that α is not 1 − 1 on the integer points of C. See [16] for several versions of Siegel’s bound. (ii) Show that the exponent n/(m − n) cannot be improved. □

Exercise 9.4: Show:

(i) If d_{1}, d_{2} are co-prime then ℤ_{d_{1}} ⊕ ℤ_{d_{2}} ≈ ℤ_{d_{1}d_{2}}.

(ii) Every finite cyclic group is a direct sum of cyclic groups of prime power.

(iii) Every finitely generated Abelian group written as a direct sum G = ⊕_{i=0}^{ℓ}H_{i} where H_{0} is a free Abelian group of rank β, and H_{i} ≈ ℤ_{q_{i}} where q_{i} is a prime power (i = 1, …, ℓ). Moreover, the numbers β, q_{1}, …, q_{ℓ} are uniquely determined by G. (These q_{i}’s are called the invariant factors of G. □

References

[1] A. Bachem and R. Kannan. Polynomial algorithms for computing the Smith and Hermite normal forms of an integer matrix. SIAM Journal Computing, 8:499–507, 1979.

[2] E. F. Bareiss. Sylvester’s identity and multistep integer-preserving Gaussian elimination. Math. Comp., 103:565–578, 1968.

[3] E. F. Bareiss. Computational solutions of matrix problems over an integral domain. J. Inst. Maths. Appl., 10:68–104, 1972.

[4] T. J. Chou and G. E. Collins. Algorithms for the solution of linear Diophantine equations. SIAM Journal Computing, 11:687–708, 1982.

[5] J. H. Davenport, Y. Siret, and E. Tournier. Computer Algebra: systems and algorithms for algebraic computation. Academic Press, 1988.

[6] Jack Edmonds. Systems of distinct representatives and linear algebra. J. Res. National Bureau of Standards, 71B:241–245, 1967.

[7] F.R. Gantmacher. Theory of Matrix, volume 1. Chelsea Publishing Co., New York, 1959.

[8] J. L. Hafner and K. S. McCurley. Asymptotically fast triangularization of matrices over rings. SIAM Journal Computing, 20:1068–1083, 1991.

[9] L. K. Hua. Introduction to Number theory. Springer-Verlag, 1982.

[10] C. S. Iliopoulous. Worst-case complexity bounds on algorithms for computing the canonical structure of finite Abelian groups and Hermite and Smith normal form of an integer matrix. SIAM Journal Computing, 18:658–669, 1989.

[11] Tudor Jebelean. An algorithm for exact division. Journal of Symbolic Computation, 15(2):169– 180, 1993.

[12] D. E. Knuth. The Art of Computer Programming: Seminumerical Algorithms, volume 2. Addison-Wesley, 2nd edition edition, 1981.

[13] Heinz Lüneburg. On the computation of the smith normal form. Preprint 117, Universitaet Kaiserslautern, Fachbereich Mathematik, Erwin-Schroedinger-Strasse, D-6750 Kaiserslautern, Germany, March 1987.

[14] Morris Newman. Integral Matrices. Pure and Applied Mathematics Series, vol. 45). Academic Press, New York, 1972.

[15] T. Sasaki and H. Murao. Efficient gaussian elimination method for symbolic determinants and linear systems. ACM Trans. on Math. Software, 8:277–289, 1982.

[16] Wolfgang M. Schmidt. Diophantine Approximations and Diophantine Equations. Lecture Notes in Mathematics, No. 1467. Springer-Verlag, 1991.