Finding exact rational solution to linear integer equations in Matlab

Question

I have a linear system of equations $$Ax=b$$ where $A$ is an $N\times N$ matrix with integer values, and $b$ is a $N\times 1$ vector with integer values. Due to prior knowledge, I know that I am guaranteed that there exists exactly one rational solution $x$ such that $x=A^{-1}b$ (as well as that $A$ has full rank).

I am now searching for an algorithm in Matlab (or any other language which I can rewrite into Matlab) which provides the exact rational solution to this problem for a given $A$ and $b$, for example in the form $\vec{x}=\frac{r}{s}\vec{t}$, where $r$ and $s$ are scalar integers and $t$ is an $N\times 1$ vector of integers, or in the form (Matlab notation) x=r./s where r and s are $N\times 1$ integer vectors.

I have tried [r,s]=rat(A\b), which works for simple enough cases but quickly results in rounding problems (i.e. r./s only approximates $x$, but I need the exact solution). Using symbolic calculations works; however, I have to compile the program to stand-alone, and the Matlab Compiler does not support the Symbolic Toolbox as far as I know and as far as I tried.

Some stats: $N$ should be allowed to be at least around $100$ or higher. The values in $A$ can become rather large, such that I am already now using int64. I understand that at one point I will run into numeric problems, and it's sufficient for me to postpone these kinds of problems as long as possible (I am grateful for any solution, even if it only works for smaller $N$). However, I need the precise solution, and I prefer an error message over any approximation, no matter how good. Runtime of the algorithm is secondary, and it's OK if it needs an hour or so for $N=100$ (being more efficient is however a plus).

Answer 1

You might look up "complete fraction-free factorization" methods. The paper "Generalized fraction-free LU factorization for singular systems with kernel extraction" contains pseudo-code.

Answer 2

We have the linear system $\rm A x = b$, where $\mathrm A \in \mathbb Z^{n \times n}$ and $\mathrm b \in \mathbb Z^n$. We would like to find rational solutions $\mathrm x \in \mathbb Q^n$. Computing the Smith normal form of matrix $\rm A$ in MATLAB:

[U, V, S] = smithForm(A) 

we obtain unimodular matrices $\rm U$ and $\rm V$ (i.e., integer matrices that have integer inverses) and diagonal matrix $\rm S = U A V$. Hence, $\rm A x = b$ can be rewritten as $\rm S V^{-1} x = U b$. Let $\rm y := V^{-1} x$. We then obtain the linear system $\rm S y = U b$, which is easy to solve because $\rm S$ is diagonal. Since both $\rm S$ and $\rm U b$ are integer-valued, the solution will be a rational vector. Lastly, we compute $\rm x = V y$, which is also a rational vector because $\rm V$ is also integer-valued.

Answer 3

The problem can be solved using the Matlab symbolic toolbox:

N = 100; %Dimensionality
A = randi(1000, [N, N]);
assert(rank(A) == N);
B = randi(1000, [N, 1]);

As = sym(A);
Bs = sym(B);

tic
xs = linsolve(As, Bs); % is a rational number
toc

disp(xs);

xs is represented as a vector of rational numbers. The runtime on an i5 is about 9 seconds for the linsolve call.