Skip to main content.

Solving Certain Linear Systems over the Integers.

.

Assume you have a linear system of equations A x = b, where A \in \Z^{n \times n} and b \in \Z^n. Assume that \det A \neq 0, and that we know that a solution in \Z^n exists. One question is: how can we efficiently compute x? Clearly, any algorithm solving linear systems over the integers or rationals will do; for example, the algorithms from the Integer Matrix Library by Z. Chen, C. Fletcher and A. Storjohann will do. That library will find any solution x \in \Q^n, and also does not require that A is invertible (over the rationals) or that A is square. But for our purposes, using such a general solver is overkill.

Note that the below material is well-known among experts.

Let p be any prime not dividing \det A. Then A modulo p is invertible, and modulo p, the system A x = b has a unique solution. Moreover, for any integer e, the system A x = b has a unique solution modulo p^e: this is true since A is also invertible modulo p^e – for that, it suffices to check that the determinant of A is a unit, which it is since it is coprime to p^e. Moreover, if y \in \Z^n is a solution to A x = b over the integers, then y modulo p^e is the unique solution of A x = b modulo p^e. Hence, if we choose e such that \frac{1}{2} p^e bounds all coefficients of the solution y, we can recover a solution to A x = b over the integers from a solution to A x = b modulo p^e, by chosing the unique preimages in (-\tfrac{1}{2} p^e, \tfrac{1}{2} p^e].

This opens the question on how to solve A x = b modulo p^e. For that, a Hensel-like lifting technique can be used. (In fact, this follows from Bourbaki's generalization since the Jacobian of the map f : (\Z/p^e\Z)^n \to (\Z/p^e\Z)^n, x \mapsto A x - b equals A.) Assume that we have an x \in \Z^n which satisfies A x \equiv b \pmod{p^{e-1}}. We want to find x' \in \Z^n with A x' \equiv b \pmod{p^e}. Write x' = x + p^{e-1} x'' with x'' \in \{ 0, \dots, p - 1 \}^n. As A x' = A x + p^{e-1} A x'', and as A x - b is divisible by p^{e-1}, we obtain the linear system A x'' \equiv \frac{A x - b}{p^{e-1}} \pmod{p}. Hence, it suffices to solve e linear systems over the prime field \F_p to solve A x = b over \Z/p^e\Z.

This yields the following algorithm:

  1. Choose p := 2.
  2. Solve A x = b modulo p.

  3. If a unique solution exists:

    1. Set e = 0 and lift x to the integers with coordinates in (\tfrac{1}{2} p, \tfrac{1}{2} p].
    2. Compute c := A*x - b.
    3. If c = 0, return x.
    4. Solve A y = c/p^e modulo p.
    5. Set x := x + y*p^e and e := e + 1.
    6. Adjust x modulo p^{e+1} such that all coefficients are in (\tfrac{1}{2} p^{e+1}, \tfrac{1}{2} p^{e+1}].
    7. Go back to Step 3.2.

    Else:

    1. Choose the next prime p and go back to Step 2.

The only subprogram we need is a linear systems solver for A x = b with square A over a finite field, which returns information on the number of solutions. (Note that \det A is not divisible by p if and only if there is a unique solution.) If more information is known on the matrix A, for example its determinant has been already computed, this information can be used as well.

Let us analyze the running time of this algorithm. Denote by NP(A) the smallest prime not dividing \det A, and by S(n, p) the time the linear system solver over \F_p needs to solve a system of size n. Let \|A\|_\infty (resp. \|x\|_\infty resp. \|b\|_\infty) denote the largest absolute value of an coefficient of A (resp. x resp. b).

Clearly, the number of iterations is in O(\log_{NP(A)} \|x\|_\infty) = O(\frac{\log \|x\|_\infty}{\log NP(A)}). In each iteration, one linear system over \F_p of size n has to be solved, and A x - b has to be evaluated. The former takes S(n, NP(A)) operations, and the latter involves n^2 multiplications and additions of integers of size O(\log \|A\|_\infty) and O(e \log NP(A)), and n substractions of integers of size O(\log \|A\|_\infty + e \log NP(A)) and O(\log \|b\|_\infty). For simplicity, assume that \log \|b\|_\infty = O(\log \|A\|_\infty). Finally, to compute x = x + y p^e, we need n multipliations of integers of size O(\log NP(A)) and O(e \log NP(A)), and n additions which can be neglected. Clearly, the n multiplications can also be neglected, since the evaluation of A x already is slower.

Let M(m) denote the time a multiplication of two numbers of size m needs. Then inside the main loop, we need

O\bigl(S(n, NP(A)) + n^2 M(\max\{ \log \|A\|_\infty, e \log NP(A) \})\bigr)

time units, and the main loop alltogether needs

& O\Biggl(\sum_{e=1}^{\frac{\log \|x\|_\infty}{\log NP(A)}} \biggl( S(n, NP(A)) + n^2 M(\max\{ \log \|A\|_\infty, e \log NP(A) \}) \biggr) \Biggr) \\ {}={} & O\Biggl(\frac{\log \|x\|_\infty}{\log NP(A)} \bigl( S(n, NP(A)) + n^2 M(\max\{ \log \|A\|_\infty, \log \|x\|_\infty \}) \biggr) \Biggr)

time units. Finding p needs

O\Biggl(\frac{NP(A)}{\log NP(A)} S(n, NP(A)) \Biggr)

time units (using the Prime Number Theorem).

Assuming that we use a naive Gaussian algorithm as well as naive multiplication, i.e. S(n, p) = n^3 (\log p)^2 and M(m) = m^2, we obtain a total running time of

O\Biggl( & n^3 \bigl( \log \|x\|_\infty + NP(A)\bigr) \log NP(A) \\ & {}+ n^2 \max\biggl\{ \frac{(\log \|A\|_\infty)^2 \log \|x\|_\infty}{\log NP(A)}, \frac{(\log \|x\|_\infty)^3}{\log NP(A)} \biggr\} \Biggr).

Using fast multiplication, i.e. M(m) = m^{1 + \varepsilon} (using FFT methods), and fast linear system solving, i.e. S(n, p) = O(n^\omega (\log p)^{1 + \varepsilon}), where \omega \le 2.376, we obtain a total running time of

O\Biggl( & (NP(A) + \log \|x\|_\infty) n^\omega (\log NP(A))^{\varepsilon} \\ & {}+ n^2 \max\biggl\{ \frac{\log \|x\|_\infty (\log \|A\|_\infty)^{1+\varepsilon}}{\log NP(A)}, \frac{(\log \|x\|_\infty)^{2 + \varepsilon}}{\log NP(A)} \biggr\} \Biggr)

Now let us try to eliminate NP(A) from this expression. Clearly, the the second part, we can use that NP(A) \ge 2. To eliminate NP(A) from the first part, we need to find an upper bound. For that, let us first stick to NP(t), the smallest prime not dividing the integer t. (Letting A be a 1 \times 1-matrix A yields NP(t) = NP(A); in general, NP(A) = NP(\det A) using this notation.) Now t is divisible by \prod_{p < NP(t)} p, whence for t < \prod_{p < x} p we have NP(t) < x. Note that for integral x, \log \bigl( \prod_{p < x} p \bigr) = \vartheta(x - 1) \le \vartheta(x) \sim x, where \vartheta denotes the Chebyshev function. Using known bounds on \vartheta(x), we get

\prod_{p < x} p = \exp(x + O(x/\log x)) = \exp((1 + o(1)) x).

Therefore, \prod_{p < x} p > \exp((1 - \varepsilon) x) becomes true for x \to \infty for every \varepsilon. This shows that NP(t) < \frac{\log t}{1 - \varepsilon} eventually holds for t \to \infty, yielding NP(t) = O(\log t) and, thus, NP(A) = O(\log \det A). Using the Leibniz formula, \log \det A = O(n \log n + n \log \|A\|_\infty).

Finally, we can use some linear algebra to bound \|x\|_\infty in terms of A and b. First note that A A^\# = (\det A) I_n, where I_n denotes the n \times n identity matrix and A^\# denotes the adjungate matrix of A. As x = A^{-1} b = \frac{1}{\det A} A^\# b, we see that it suffices to bound \|A^\#\|_\infty. Now the coefficients of A^\# are determinants of (n - 1) \times (n - 1) matrices with coefficients coming from A, whence \|A^\#\|_\infty \le (n - 1)! \|A\|_\infty^n. Therefore,

\log \|x\|_\infty \le n \log n + n \log \|A\|_\infty + \log \|b\|_\infty = O(n \log \|A\|_\infty)

when assuming that \log n, \log \|b\|_\infty = O(\log \|A\|_\infty).

This can be combined into the following theorem:

Theorem.

Assuming that \log n = O(\log \|A\|_\infty) and \log \|b\|_\infty = O(\log \|A\|_\infty), the above algorithm needs

O\bigl( n^5 (\log \|A\|_\infty)^3 \bigr)

time units to compute the unique solution of A x = b using naive arithmetic in \Z, \F_p, and naive Gaussian elimination to solve linear systems over \F_p. Using fast linear algebra and fast multiplication, we only need

O\bigl( n^{4 + \varepsilon} (\log \|A\|_\infty)^{2 + \varepsilon} \bigr)

time units for any \varepsilon > 0.

Categories: Computational Number Theory Linear Algebra
Tags:

Comments.

No comments.