From: rusin@eiger.math.niu.edu (Dave Rusin)
Newsgroups: sci.math
Subject: Re: Euclidean Algorithm For Factoring Polynomial Matrices... How?
Date: 1 Feb 1997 07:49:11 GMT

In article <5c9jnl$2f0e@uvaix2e1.comp.UVic.CA>,
Michael Adams <mdadams@uvaix2e1.comp.UVic.CA> wrote:
>Hi, Folks.
>
>For some DSP work that I'm doing, I need to factor an N by N matrix of
>Laurent polynomials (with N >= 3).  Each matrix factor needs to have
>the form of an elementary row/column operation.  As I understand it,
>the Euclidean algorithm can be used to do this, but I'm not quite sure
>how.  I've seen the Euclidean algorithm used to find the GCD of several
>polynomials, but I still can't quite make the connection between the
>Euclidean algorithm and the factorization problem I mentioned above.
>
>Can anyone provide any pointers?  A brief explanation, some pseudocode,
>a worked-through example, or real code would certainly help a lot.

You need to build on both the GCD calculations and Gaussian elimination.

The connection with the latter is discussed in typical elementary
linear algebra texts: if you you pass from the matrix  A  to a
row-echelon form  B (for example) by a sequence of elementary
row operations, then apply each of these row operations (separately) 
to the identity matrix to obtain elementary matrices  E_i. Then
	B = E_n * ... * E_2 * E_1 * A.
Since the  E_i are invertible (and the inverses also elementary) we
have a factorization  A = E_1^(-1) * E_2^(-1) * ... * E_n^(-1) * B.
If you use column operations as well, then elementary matrices will
appear to the right of  B  too. If  A  is a nonsingular square
matrix of real numbers,  B  can be taken to be the identity 
and so need not appear. I assume the resulting product is
the kind of factorization you seek.

Now imagine you wish to mimic just this procedure on your matrix of
Laurent polynomials. Let's see, you begin in the upper left; interchange
rows if necessary to make that entry nonzero. Then divide through by that
entry to get a 1, then subtract off multiples to reduce the rest of that
column to zeros. But wait, in the ring of Laurent polynomials, not every
nonzero element is invertible. Thus you must accept the fact that
the pivot entry will be some nonzero element of the ring. You _can_
divide through to make sure that upper-left entry is a monic polynomial,
and indeed you may assume all entries in the first column are, too.

But now how can we arrange to nearly clear the first column? Suppose
for example the first column's entries are   x^2+3x+2, x^2+2x+1, and zero.
You know the routine for computing the GCD of the first two polynomials:
subtract and so replace the second by x+1, then subtract a multiple of
this from the first to get 0; arriving at  0  informs us that the last
term (x+1) is the GCD.

Well, since you're just subtracting multiples, you can just as well do
this to all entries in a row, which is an ordinary row operation. In 
this way you see that by doing row operations you can arrange it so that
the one nonzero entry in the first column is just the  GCD of the terms
which were there in the first place. In particular, you'll have a
matrix whose upper-left entry is now of no greater degree than before.
('degree' for Laurent polynomials is the difference between the greatest
and least exponents used).

Now, what you do next is perhaps not quite what you'd expect. You see,
continuing in this way you can write  A  as a product of elementary 
matrices times  B, the matrix (in row-echelon form) which results from
performing row-reduction. This  B  will be upper-triangular. Here
you haven't quite got the result you've asked for, but then again you've
only used row operations. If you really need  B  to be a diagonal
(obviously then a product of elementary matrices) you'll need to do
column operations too. But don't wait until all the row operations are done.

Instead, after you've cleared the first column, use column operations
to clear the first row, replacing the upper-left entry by the GCD of
the entries in the first row. This must divide what was there before,
and so be of smaller degree -- unless it's of the same degree, in
which case the old upper-left entry _is_ the GCD (in that case,
subtracting simple multiples of the first column from the others will
clear the first row).

Iterating this procedure until the GCD stops decreasing, you finally
end up with a matrix which has a small entry in the upper left, and
zeros in the rest of the first row and first column. Thereafter, you
can use the same algorithm on the resulting (N-1) x (N-1) matrix to
the bottom right.

Continue in this way until a diagonal matrix (clearly a product of
elementary matrices) results.

Incidentally, we never use any row operations which change the
determinant by more than a unit, so the determinant of the original
matrix is (up to a unit) the product of the terms on the final diagonal.
Moreover, one can show that the diagonal terms obtained are the GCDs
of the whole matrix from which they were obtained, and so divide
each other as you progress down the diagonal. They are invariants
of the original matrix (up to units).

You'll find discussion of this in any nice abstract algebra book
(look for "modules over a PID", "finitely-generated abelian group structure"
or "Jordan canonical form").

As a practical matter, you want to use row- and column- operations more
or less simultaneously. For example, one usually begins by rotating into
the upper left entry that entry from the first row which is of lowest
degree; well, why not go whole hog and use row- _and_ column- swaps to
put into that position the one entry from the whole matrix which is of
lowest degree. (Interestingly, the final results can be achieved by
doing only row operations first, then only column operations -- this is
a consequence of associativity of matrix multiplication -- but you'll
never _deduce_ optimal sets of operations in this way.)

dave