Week of 8/4 Update

This week I've been working on on fleshing out the algorithm I described last week for computing the last invariant factor of a polynomial matrix. I've written a high-level implementation in Sage to demonstrate proof of concept. There are a few subtleties I've uncovered:

• The degree of $\bar{x}$, $N$, needs to be bigger than I originally thought. A safe bound is $N\geq{}2n*deg(A)+deg(B)+1$ although I'm not yet sure whether this bound is sharp.
• The degree of $\vec{b}$, $d_{b}$, can be fairly small, about $2\log_{q}n$. This ensures that the components are probably distinct.
• When computing $\frac{n_{i}}{d}$ such that $\frac{n_{i}}{d}\equiv{}\bar{x_{i}}(mod\;x^{N})$ there's one such for any given choice of $k=deg(d)$. However only one such choice of $k$ will actually solve the system. To figure out which choice is correct I run the extended euclidean algorithm for several rows of $\bar{x}$ simultaneously and mark values of $d$ that occur in more than one row as candidates.
• Given a candidate $d$ the other $n_{i}$'s can be determined quickly by dividing $\frac{\bar{x}_{i}}{d}$
• Any particular $\bar{x}_{i}$ has a chance of producing the wrong $d$ due to cancellation in the Cramer's rule formula $\bar{x}_{i}=\frac{|A_{i}|}{|A|}$. This occurs with probability $\frac{1}{q}$. Checking more $\bar{x}_{i}$'s for candidate $d$'s fixes this. To ensure that we get the correct $d$ as a candidate with probability $p$ we need to check $qbinom(p,n,1-\frac{1}{q})$ where qbinom is the binomial quantile function.
• Sometimes there will be no $\alpha$ such that $|A(-\alpha)|\neq{}0$, even if $A$ is non-singular (over the rational functions). In this case we can still find an answer by finding a random solution to the Block-Toeplitz system: $\begin{pmatrix}A_{d_{b}}&A_{d_{b-1}}&...&A_{0}&0&...\\A_{d_{b+1}}&A_{d_{b}}&...&A_{0}&0&...\\&&\ddots{}&&\\0&...&A_{d_{A}}&A_{d_{A}-1}&...&A_{0}\end{pmatrix}\begin{pmatrix}\bar{x}^{(0)}\\\bar{x}^{(1)}\\\vdots{}\\\bar{x}^{N}\end{pmatrix}=\vec{0}$

 Next week I'll start wrapping things up as the Google Summer of Code draws to a close.

Week of 7/28 Update

This week I've started looking into an alternative way of computing the modulus for the polynomial ring used by SmithFormIliopolus. The problem with PolyDet is that it's painfully slow, mostly due to its heavy use of arithmetic over field extensions. The new approach would avoid the use of extensions and is based on a paper by Dixon (Exact Solutions of Linear Equations Using P-Adic Expansions (1982)). The difference is that here we're working with polynomial, not p-adic representations of integers. The goal is to compute the largest invariant factor of a matrix A. For our purposes, this is actually better than the determinant since it keeps our intermediate polynomials smaller. Here's the algorithm:

• Choose a random polynomial vector $\vec{b}$ over the ground field $\mathbb{F}$
• The degree of $\vec{b}$ (defined as the maximum of the degrees of the elements of $\vec{b}$) should be large enough to sample the range of $A$, but how big exactly remains unclear. 
• Choose $\alpha$ to make $A(-\alpha)$ non-singular then let $A\leftarrow{}A(x-\alpha)$
• Solve $A\vec{x}=\vec{b}$ for the first $2n$ terms of a power series solution for $\vec{x}$ for $n$ the sum of the degrees of each row of $A$
• Perform a rational reconstruction to solve $\vec{x}=\frac{\vec{n}}{d}$ for $d$
• Take a few random linear combinations of the elements of $\vec{x}$
• For each linear combination find $f$ and $g$ of minimal degree such that $g\vec{x}\equiv{}f$ (mod $x^{2n}$)
• Take the least common multiple of the various $g$'s
• Return the polynomial $d$

Asymptotically this algorithm require about the same number of field operations as PolyDet. However, because it only works over the ground field I'm hoping it will be much faster. I'm working on a Sage example of this algorithm now to work out the details. Hopefully I'll then be able to get a Linbox implementation up and working to compare against PolyDet.