5
$\begingroup$

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

|cite|improve this question
$\endgroup$
  • $\begingroup$ One thing to try would be implementing Gaussian elimination symbolically yourself. $\endgroup$ – tch Feb 16 at 15:55
  • $\begingroup$ Have you tried the Smith or Hermite normal forms? $\endgroup$ – Rodrigo de Azevedo Feb 17 at 13:12
  • $\begingroup$ @RodrigodeAzevedo What exactly should I do with these forms? Respectively, how to I get to rational solutions once I have these forms? $\endgroup$ – Engineer trying math Feb 20 at 10:58
  • $\begingroup$ Sounds like you're getting to the point at which Matlab becomes more a hindrance than a useful tool. $\endgroup$ – Federico Poloni Jul 2 at 18:10
3
$\begingroup$

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.

|cite|improve this answer
$\endgroup$
1
$\begingroup$

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.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ Thanks a lot. So, the good news is that I implemented this algorithm and it works. The bad news is that the function smithForm needs the symbolic toolbox and is thus not compatible with the Matlab compiler (and there are easier ways to solve the problem as soon as we can use the symbolic toolbox). I tried to replace it by some other implementation of the smith normal form, but most of them I tested seem to be either buggy or to have problems handling big numbers. I currently try to find something which I can make stable enough to use, let's see... $\endgroup$ – Engineer trying math Feb 24 at 19:47
  • $\begingroup$ @Engineertryingmath Yes, big numbers. I am not surprised. Take a look at this. $\endgroup$ – Rodrigo de Azevedo Feb 25 at 9:36
  • $\begingroup$ @Engineertryingmath This is also interesting. $\endgroup$ – Rodrigo de Azevedo Feb 25 at 12:30
  • $\begingroup$ The main problem is not computational time, but that the values during the computation quickly approach flintmax('double') (beyond which $\endgroup$ – Engineer trying math Feb 26 at 17:20
  • $\begingroup$ The problem is not computational time, but that the values during the computation of the smith normal form for the algorithms I am looking at quickly approach flintmax('double') (beyond which we get numerical errors). So, I rewrote one of these algorithms to use variable precision integers instead of doubles (using the vpi toolbox). However, while this is working in principle, everything gets incredibly slow. The profiler says that this is due to matrix multiplications, which seem to be implemented rather inefficiently for VPIs. Solving one problem just seems to create the next one... $\endgroup$ – Engineer trying math Feb 26 at 17:29
0
$\begingroup$

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.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ Thanks a lot. However, as I wrote in my question, I cannot use symbolic calculations, since the Symbolic Toolbox is not supported by the Matlab Compiler. Lines 6 and 7 thus will return errors, since the function "sym" does not exist in deployed code. $\endgroup$ – Engineer trying math Feb 16 at 12:37
  • $\begingroup$ Sorry, I overlooked this Info. You could maybe use SymPy and distribute it with your program (the tarball is smaller than 10 MB). However this relies on having python on the target platform. You could maybe also use GMP (gmplib.org) and implement a basic linear system solver using the rational numbers provided there. You could then compile and link to matlab. $\endgroup$ – Lukas Koestler Feb 16 at 12:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.