computing the inverse

Why

For many calculation algorithms (division of "BigInts", division of polynomials, division of power series, ...) it is desirable to compute efficiently the (multiplicative) inverse of an element b of some ring R. Here's a fairly general method, that works in many real problems. Like the Babylonian square root algorithm, it can be derived from the Newton-Raphson method; also like it, you don't have to use the iteration.

How

We assume we have some norm |.|:R→Q⁺ (or sometimes to the real numbers), i.e. just a function satisfying

|1|=1
|p*q|≤|p|*|q|

Now, if you've glanced at the Newton-Raphson method, you'll know that it requires division. How can we use an algorithm which requires division in order to implement division? We'll arrange matters so all divisions in the algorithm cancel out, leaving us with just multiplications and additions.

Let b∈R; we shall compute 1/b. We decide to solve the equation 1/x-b=0. Dumbly applying the rules for deriving the Newton-Raphson iteration (after all, derivatives might not mean anything in our ring R!), we see:

f(x)=1/x-b
f'(x)=-1/x²

f(x)/f'(x)=bx²-x

So our desired iteration will be

x_n+1 = x_n - f(x_n)/f'(x_n) = 2x_n-bx_n².

Why does it work?

We shall prove that this iteration converges (and rapidly, at a quadratic pace!). So we can ignore the fact that we used the Newton-Raphson method (with somewhat shaky justification) to derive our iteration -- we shall just prove ``it works''.

When can we guarantee convergence?

|x_n+1-1/b| = |2x_n-bx_n²-1/b| = |b*(x_n-1/b)²| ≤ |b|*|x_n-1/b|²

So we're guaranteed quadratic convergence if |x₀-1/b|<|b|.

Applications

Power series

Let R=Q((Z)) be the Laurent field of formal power series with rational coefficients. A good choice of "norm" is

v(p(Z)) = the lowest power of Z appearing in p(Z)
|p(Z)| = 2^-v(p(Z))

"Convergence" in this norm means agreement of increasing numbers of low-order terms.

We show how to compute the inverse of a power series of norm 1. For general p(Z), just note that we can express the inverse based on the inverse of such a power series:

1/p(Z) = Z^v(p(Z))*1/(Z^-v(p(Z))p(Z)).

And if v(p(Z))=1, it turns out that a good initial guess is p(0), the coefficient of the 1's term in p(Z).

Let us compute the inverse of Z²+1. The first few iterations:

1. a₀(Z) = 1 (the coefficient of the 1's term).
2. a₁(Z) = 2a₀(Z)-(Z²+1)a₀(Z)² = 1-Z²
3. a₂(Z) = 2a₁(Z)-(Z²+1)a₁(Z)² = 1-Z²+Z⁴-Z⁶

Indeed, every iteration doubles the number of correct terms! (The correct answer is 1-Z²+Z⁴-Z⁶+Z⁸-...).

Polynomials mod Z^N

A more practical setting is used when dividing polynomials. We want to find the multiplicative inverse of a polynomial b(Z) modulo Z^N (for more information, see the full algorithm under polynomial division). This exists (when working over a field; for polynomials over a ring we require b(0) invertible...) iff b(0)≠0.

We proceed exactly as in the power series case. Except that now there cannot be norms smaller than 2^-N. So we may conclude our iteration after log₂(N) iterations (rounded up), knowing we hold the precise answer.

Integers mod 2^N

Suppose you're trying to divide huuuge numbers. Computer arithmetic is usually performed modulo 2^N, so we shall do so here, too. The 2-adic norm will be our norm; we can invert any odd number.

For instance, to invert 13 modulo 2^16:

1. a₀=1 (≥1 bit of accuracy)
2. a₁=-11≡65525 (mod 2^16)=1111111111110101₂ (≥2 bits of accuracy; in fact, 4 bits!)
3. a₂=-1595≡63941 (mod 2^16)=1111100111000101₂
4. a₃≡20165=100111011000101₂ is the correct answer.

Matrices

Some square matrices have an inverse. We can compute it rapidly. Suppose we want to invert a matrix of real or complex numbers B. If the matrix norm |I-B^-1|<1, then by the argument above the algorithm, started from A₀=I, will converge.

But many ("most"!) nonsingular matrices B are not of this form. Still, this is the basis for a very strong method of inverting a matrix.