Suppose that I have a linear system $AX=b$ with $A\in\mathbb{Z}^{n\times m}$ and $b\in\mathbb{Z}^{n}$. Assume that $AX=b$ has exactly one rational solution. Then how can I obtain this solution efficiently?


I know LU decomposition is a choice for this purpose. But it seems that the denominators may expand in the process of computing a LU decomposition, which I hope to avoid if the rational solution has small denominators.


If I already have a numerical solution of $AX=b$, does it help?



How to obtain the rational solution of a linear system efficiently?
The following answer is only theoretical. Suppose you first reduce your matrix to [Smith normal form](https://en.wikipedia.org/wiki/Smith_normal_form) $D = S A T$, where $S$ and $T$ are integral matrices with integral inverses (the determinants of $S$ and $T$ are each either $+1$ or $-1$) and $D$ is diagonal, whose diagonal elements are the so-called *elementary divisors* of $A$. Then you reduce your problem to solving $DY=c$, where $c = S b$ and $X = T Y$. Since by $c$ is still integral and $D$ is now diagonal, the denominators in $Y$ will be no larger than the largest entry in $D$ (the largest elementary divisor of $A$). Then, since $X$ is obtained from $Y$ by multiplication with another integral matrix, the denominators in $X$ will be no larger than the lcm of the denominators in $Y$.


Unfortunately, I don't think that reduction to Smith normal form is a cheap operation, so whether it is worth going through that step in your problem will depend on other factors.


EDIT: Actually, the Smith normal form is overkill for your problem. The [Hermite normal form](https://en.wikipedia.org/wiki/Hermite_normal_form) is basically the integer version of the LU decomposition that you already had in mind: $U = S A$, where $U$ is integral upper triangular, and $S$ is integral with integral inverse. The new system you need to solve is $U X = S b$, which can be done by back substitution, where at each step you only divide by the corresponding diagonal element of $U$. Hence, the final denominators in $X$ will be no larger than the lcm of the diagonal elements of $U$.


The linked Wikipedia article on the Hermite normal form mentions that modern algorithms for its computation have polynomially bounded time and space complexity (respectively, in terms of the matrix size and logarithmic size of the matrix entries). If you (or your computer algebra system) use such an algorithm, it might be your best bet for *provably* avoiding intermediate expression swell. The *practicality* of this method for your problem of course would still be subject to experimentation.


EDIT 2: I've also found a very well written survey of several algorithms for solving rational linear systems, with practical timing comparisons in the case of sparse systems.



> 
> *Cook, William; Steffy, Daniel E.*, [**Solving very sparse rational systems of equations**](http://dx.doi.org/10.1145/1916461.1916463), ACM Trans. Math. Softw. 37, No. 4, Paper No. 39, 21 p. (2011). [ZBL1365.65121](https://zbmath.org/?q=an:1365.65121). (author's [copy](https://www.math.uwaterloo.ca/~bico/papers/rational.pdf))
> 
> 
> 


