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 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