MATHEMATICS OF COMPUTATION VOLUME 000, NUMBER 0 0000, PAGES 000–000

MATRIX INVERSION MADE DIFFICULT

DANIEL J. BERNSTEIN

Abstract. If ADA = A and DAD = D then D is a quasiinverse of A. This paper presents generalized Gaussian elimination, which builds a quasiinverse of an arbitrary homomorphism out of quasiinverses of two smaller homomorphisms. In conjunction with fast matrix multiplication, generalized Gaussian elimination is a fast, reasonably stable method that can (1) invert any invertible matrix, (2) solve any soluble set of linear equations, and (3) compute a basis for the kernel of any matrix. A fast algorithm for (1) was already obtained by Bunch and Hopcroft, but the construction here is simpler. Fast algorithms for (2) and (3) appear to be new.

1. Introduction

Let A : V → W and D : W → V be module homomorphisms. Then D is a quasiinverse of A if ADA = A and DAD = D.

If A has a left inverse then its quasiinverses are exactly its left inverses, and if A has a right inverse then its quasiinverses are exactly its right inverses. In particular, when A is invertible, its only quasiinverse is its inverse.

Any quasiinverse D determines a coimage of A in V, namely im D, and a cokernel of A in W, namely ker D. (It is tempting to say that D is co A; D reverses all exact sequences related to A.) Then V splits as (ker A) ⊕ (im D) and W splits as (ker D) ⊕ (im A). Note that A is an isomorphism from im D to im A, with D as its inverse.

A quasiinverse D respects a basis for W if ker D is generated by a subset of that basis. Similarly, D respects a basis for V if im D is generated by a subset of that basis. Note that the corresponding columns of A then form a basis for im A.

Any matrix over a field has a quasiinverse respecting any given bases for V and W.

Outline of this paper. Section 2: Generalized Gaussian elimination constructs a quasiinverse of A given quasiinverses for a pivot and a reduced map.

Section 3: Given a matrix A over a field, Algorithm Q computes a respectful quasiinverse D of A, along with bases for ker D and im D. Algorithm N computes a basis for ker A. Algorithm S solves linear equations.

Section 4: The algorithms of section 3 run at essentially the speed of matrix multiplication.

\begin{matrix}c \\ ◯ \end{matrix}0000 American Mathematical Society 0025-5718/00 $1.00 + $.25 per page

Received by the editor 19951021 (draft 7).

1991 Mathematics Subject Classification. Primary 65F05, 15A09.

Thanks in advance to anyone who suffers through the first few drafts of this paper.

Fast matrix operations. Strassen [4] showed that one can multiply two n-by-n matrices in time n^{β} for a certain β < 3. (Currently the best known β is a bit below 2.38 [2].) For moderately large matrices Strassen’s method is much faster than the usual algorithm.

Strassen then suggested block Gaussian elimination as a fast method to invert an n-by-n matrix. Unfortunately, this method can fail even if the original matrix is invertible. (It is also rather unstable.) Bunch and Hopcroft [1] fixed this problem with a tricky series of nested pivots; they proved that one can invert any invertible n-by-n matrix at roughly multiplication speed.

Generalized Gaussian elimination achieves the same result as the Bunch-Hopcroft method for invertible matrices, in a much less convoluted way. Strassen, Bunch, and Hopcroft chopped up A both horizontally and vertically; generalized Gaussian elimination chops up A in just one direction. This simplification depends crucially on the concept of a quasiinverse.

Generalized Gaussian elimination also produces useful results for non-invertible matrices. Examples: (1) it computes a basis for ker A; (2) it points out a basis for im A; (3) it solves systems of linear equations—if D is any quasiinverse of A, and w ∈ im A, then w = A(Dw).

Relation to Gaussian elimination. In section 2, S may be chosen to be onedimensional. This special case is exactly Gaussian elimination with partial pivoting [XXX]. The latter algorithm is normally viewed as searching for a pivot element in (let’s say) the first row of A, and then using the corresponding column to reduce A. I^{1} prefer to view the entire first row as being the pivot. The inverse of the pivot element is really a quasiinverse of the pivot row. In generalized Gaussian elimination, the pivot may be any subset of the rows of A; the quasiinverse of the pivot amounts to the inverse of some maximal nonsingular square submatrix.

Quasiinverses in maximum generality. I have stated the results of section 2 for modules over any ring, not just vector spaces over a field, on the theory that unnecessary hypotheses distract attention from the essence of the proofs. For a more thorough application of this philosophy, see [3].

One could go even further, and rephrase section 2 in the language of abelian category theory; note that diagram theorems true for modules are true for all abelian categories [XXX]. I have no intention of carrying out this project.

Quasiinverses in the literature. It seems that there is no standard name for my concept of a quasiinverse. The word “quasiinverse” has been used [XXX] for a map D satisfying merely ADA = A. I suggest that such a map be called a “right quasiinverse” or “splitting” of A.

Any matrix over the real numbers has a unique pseudoinverse, a particular quasiinverse determined by the singular value decomposition. The singular value decomposition supplies orthonormal bases for im D and im A upon which A acts diagonally. See [XXX] for further details.

Another example of a quasiinverse is [XXX]; thanks to Henri Gillet and Lawrence Ein for bringing this to my attention.

Note that the construction of the reduced map C from the pivot B in Theorem 2 generalizes the construction of the Schur complement [XXX].

1I have abandoned the practice of mechanically replacing “I” with “we” in my papers.

Terminology. I abbreviate module homomorphism as map. Following convention and ignoring psychology I write maps on the left. S ⊕ T is the direct sum of S and T; it is the same as S + T, with the side condition that S ∩ T = {0}. #Z is the number of elements of Z.

2. Generalized Gaussian elimination

Consider a map A : V → S ⊕ T. This section presents generalized Gaussian elimination, which builds a quasiinverse D of A out of (1) a quasiinverse E of a pivot B : V → S and (2) a quasiinverse F of a reduced map C : V → T.

Theorem 1. Let A : V → W and π : W → S be maps. Define B : V → S by B = πA. Assume that B has a quasiinverse E : S → V. Define P : V → V by P = 1 − EB and Q : W → W by Q = 1 − AEπ. Then P E = 0, BP = 0, P P = P, QA = AP, and πQA = 0.

Proof. First P E = E − EBE = 0; next BP = B − BEB = 0 so P P = P − EBP = P; next QA = A − AEπA = A − AEB = AP; finally πQA = πAP = BP = 0. □

Theorem 2. In the situation of Theorem 1, assume that S is a submodule of W and that ππ = π. Define T = ker π. Define C : V → T by C = QA. Assume that C has a quasiinverse F : T → V. Define D : W → V by D = Eπ + P F (1 − π)Q. Then D is a quasiinverse of A. Furthermore im D = im E + im F and ker D = ker E + ker F.

The assumptions on S and W mean that W = S ⊕ T; π projects W onto S and 1 − π projects W onto T.

Proof. By Theorem 1, CE = AP E = 0 and CP = AP P = AP = C.

Now DA = EπA + P F (1 − π)QA = EB + P F C so ADA = AEB + AP F C = AEB + CF C = AEB + C = A.

Furthermore BD = BEπ + BP F (1 − π)Q = BEπ since BP = 0; thus EBD = EBEπ = Eπ. Also CD = CEπ + CP F (1 − π)Q = CF (1 − π)Q since CE = 0 and CP = C; thus F CD = F CF (1 − π)Q = F (1 − π)Q.

Hence DAD = EBD + P F CD = Eπ + P F (1 − π)Q = D.

Coimage: First DAE = EBE + P F CE = E since CE = 0; thus im E ⊆ im D. Similarly DAF = EBF +P F CF = EBF +P F = F; thus im F ⊆ im D. Conversely D = Eπ + P F (1 − π)Q = Eπ + F (1 − π)Q − EBF (1 − π)Q so im D ⊆ im E + im F.

Cokernel: If s = πs and Es = 0 then Eπs = 0 so Qs = s so Ds = Eπs + P F (1 − π)Qs = P F (1 − π)s = 0; thus ker E ⊆ ker D. If πt = 0 and F t = 0 then Qt = t so Dt = Eπt + P F (1 − π)Qt = P F (1 − π)t = P F t = 0; thus ker F ⊆ ker D. Conversely, say Dw = 0. Then Eπw = EBDw = 0, so Qw = w, so F (1−π)w = F (1−π)Qw = F CDw = 0. Thus w = πw+(1−π)w ∈ ker E +ker F.□

Theorem 3. In the situation of Theorem 2, let M be a submodule of ker A. Assume that ker C = (M + im E) ⊕ N for some module N. Then ker A = M ⊕ P N, and P is an isomorphism between N and P N.

Proof. If v ∈ N and P v ∈ M then v = P v + EBv ∈ M + im E so v = 0. Thus (1) P is injective on N and (2) M ∩ P N = {0}.

If v ∈ N then AP v = Cv = 0 so P v ∈ ker A. Hence P N ⊆ ker A.

If Av = 0 then Cv = QAv = 0 so v = m + Es + n for some m ∈ M, s ∈ S, and n ∈ N. Next A(Es + n) = Av − Am = 0, so Es + n = P (Es + n) = P n, so v = m + P n ∈ M + P N. Hence ker A ⊆ M + P N. □

3. Computation

The following algorithms use generalized Gaussian elimination, recursively, to solve various problems in linear algebra over a field.

A vector space of dimension n is normally represented inside a computer by the integer n. Its basis vectors are represented by the integers 1, 2, …, n (or 0, 1, …, n − 1). A respectful subspace is represented as a list of basis vectors. Subspace sum and intersection then correspond to list union and intersection.

Say X is a basis for V, and Y is a basis for W. Recall that a matrix is a function a that assigns a field element to each pair (y, x) ∈ Y × X. This matrix represents the map A : V → W defined by A(x) = ∑_{y} a(y, x)y. If I ⊆ X and J ⊆ Y, the notation A on J × I means the submatrix a on J × I, i.e., the function a restricted to J × I. Note that a on J × I determines A, provided that (1) the span of Y − J is contained in ker A and (2) im A is contained in the span of I.

In these algorithms, V and W are respectful finite-dimensional subspaces of some (irrelevant) larger space. A is a map from V to W. M is a respectful subspace of V, with M ⊆ ker A. X is a basis for V minus a basis for M, Y is a basis for W, and a = A on Y × X; note that a determines A. X and Y are represented as finite lists of integers.

Computing a quasiinverse. The following algorithm computes a quasiinverse for A. It avoids storing matrix rows and columns that are guaranteed by construction to be entirely zero.

Algorithm Q. Given (X, Y, a), compute Q(X, Y, a) = (I, J, D on I × J). Here D : W → V is a quasiinverse of A; I is a basis for im D, with I ⊆ X; and Y − J is a basis for ker D, with J ⊆ Y.

1. If a = 0: Set d to the empty matrix. Return ({} , {} , d).

2. If Y has just one element y: Select v ∈ X such that a(v, y) ≠ 0. Set d to the 1-by-1 matrix 1/a(v, y). Return ({v} , {y} , d).

3. Select Z ⊂ Y, Z ≠ {}. (In the language of section 2: Z is a basis for S, and Y − Z is a basis for T.)

4. Extract b = a on Z × X. (b represents B = πA.)

5. Recursively compute (I, J, e) = Q(X, Z, b). (e represents E.)

6. Extract t_{1} = a on (Y − Z) × I.

7. Compute the product t_{2} = t_{1}e, a matrix on (Y − Z) × J. (t_{2} represents (1 − π)AE.)

8. Extract t_{3} = b on J × (X − I). (t_{3} represents B off im E.)

9. Compute the product t_{4} = t_{2}t_{3}, a matrix on (Y −Z)×(X −I). (t_{4} represents (1 − π)AEB off im E.)

10. Extract t_{5} = a on (Y − Z) × (X − I). (t_{5} represents (1 − π)A off im E.)

11. Compute c = t_{5} − t_{4}, a matrix on (Y − Z) × (X − I). (c represents C = (1 − π)A − (1 − π)AEB.)

12. Recursively compute (G, H, f) = Q(X − I, Y − Z, c). (f represents F.)

13. Extract t_{6} = t_{3} on J × G.

14. Compute the product t_{7} = t_{6}f, a matrix on J × H. (t_{7} represents BF.)

15. Compute the product t_{8} = et_{7}, a matrix on I × H.

16. Compute t_{9} = t_{8} − f, a matrix on (G ∪ I) × H. (t_{9} represents −P F = EBF − F.)

17. Extract t_{10} = t_{2} on H × J. (t_{10} represents the projection of (1 − π)AE to the

coimage of F.)

18. Compute the product t_{11} = t_{9}t_{10}, a matrix on (G ∪ I) × J. (t_{11} represents −P F (1 − π)AE.)

19. Compute d = e + t_{11} − t_{9}, a matrix on (G ∪ I) × (H ∪ J). (d represents D = Eπ − P F (1 − π)AEπ + P F (1 − π).)

20. Return (G ∪ I, H ∪ J, d).

Note that the quasiinverse is square: if r is the rank of A, then #I = #J = r. Note also that it is easy to adapt Algorithm Q to compute the determinant of a maximal nonsingular submatrix of A.

Solving a system of linear equations. Fix w ∈ W. The following algorithm finds a solution z to the equation Az = w, if any solution exists. It is a straightforward optimization of the following procedure: first find a quasiinverse D of A by Algorithm Q; then compute z = Dw; then throw away D.

Algorithm S. Given (X, Y, a, w), compute S(X, Y, a) = (I, Dw on I), where D is some respectful quasiinverse of A and I is a basis for im D.

1. If a = 0: Halt if w ≠ 0; there is no solution. Otherwise set z to the empty vector and return ({} , z).

2. If Y has just one element y: Select v ∈ X such that a(v, y) ≠ 0. Set z to the length-1 vector (u/a(v, y)), where w = uy. Return ({v} , z).

3. Select Z ⊂ Y, Z ≠ {}.

4. Extract b = a on Z × X.

5. Compute (I, J, e) = Q(X, Z, b).

6. Compute t_{1}, t_{2}, t_{3}, t_{4}, t_{5}, c as in Algorithm Q.

7. Compute w_{2} = w on J.

8. Compute w_{3} = w on (Y − Z).

9. Compute w_{4} = t_{2}w_{2}, a vector on Y − Z.

10. Recursively compute (G, z_{1}) = S(X − I, Y − Z, c, w_{3} − w_{4}).

11. Extract t_{12} = b on J × G.

12. Compute w_{5} = t_{12}z_{1}, a vector on J.

13. Compute z_{2} = e(w_{2} − w_{5}), a vector on I.

14. Compute z = z_{1} + z_{2}, a vector on G ∪ I.

15. Return (G ∪ I, z).

Computing a basis for the nullspace. By assumption ker A contains a respectful subspace M of X. The following algorithm computes a (generally nonrespectful) subspace N of X such that ker A = M ⊕ N. It represents N as a matrix whose columns form a basis for N. Presumably we know a basis for M—for example, we could start with M = {0}—so we can compute a basis for ker A.

Algorithm N. Given (X, Y, a, w), compute N(X, Y, a).

1. If a = 0: Return the identity matrix on X × X.

2. If Y has just one element y: Select v ∈ X such that a(v, y) ≠ 0. For each x ∈ X with x ≠ v, form the vector x − (a(x, y)/a(v, y))v. Return the matrix of these vectors.

3. Select Z ⊂ Y, Z ≠ {}.

4. Extract b = a on Z × X.

5. Compute (I, J, e) = Q(X, Z, b).

6. Compute t_{1}, t_{2}, t_{3}, t_{4}, t_{5}, c as in Algorithm Q.

7. Recursively compute N_{1} = N(X − I, Y − Z, c).

8. Compute the product N_{2} = t_{3}N_{1}.

9. Compute the product N_{3} = eN_{2}.

10. Return N_{1} − N_{3}.

Numerical analysis. The following observations are well known for the special case of Gaussian elimination with partial pivoting.

Say A is a matrix over the real numbers. Each of the above algorithms appears to be numerically stable, provided that, in step 2, one selects the element largest in absolute value.

Inexact computations may lead to spurious failures when Algorithm S is applied to a singular system of linear equations. To avoid this one should treat sufficiently small numbers as zero in step 1. Another possibility is to skip the failure case entirely; w is thus projected (in some direction) onto im A.

4. Speed

In this section I’ll consider the speed of Algorithms Q, S, and N, from three points of view: (1) asymptotically, if we use fast matrix multiplication; (2) in the special case #Z = 1; (3) in practice.

Preliminaries. Write n = #X and m = #Y, so that a is an m-by-n matrix. In step 2 of each algorithm we split Y into two subsets, Z and Y − Z; write p = #Z, with 0 < p < m. Define r and s as the ranks of A and B respectively.

Algorithm Q computes (for m ≥ 2) five matrix products, of sizes m − p by s by s, m − p by s by n − s, s by r − s by r − s, s by s by r − s, and r by r − s by s. It also calls itself recursively with a p-by-n matrix of rank s and an (m − p)-by-(n − s) matrix of rank r − s.

Algorithm S computes matrix products of sizes m − p by s by s, m − p by s by n − s, m − p by s by 1, s by r − s by 1, and s by s by 1. It calls Algorithm Q with a p-by-n matrix of rank s; it calls itself recursively with an (m − p)-by-(n − s) matrix of rank r − s.

Algorithm N computes matrix products of sizes m − p by s by s, m − p by s by n − s, s by n − s by n − r, and s by s by n − r. It calls Algorithm Q with a p-by-n matrix of rank s; it calls itself recursively with an (m − p)-by-(n − s) matrix of rank r − s.

Fast quasiinversion. Say we can compute the product of two k-by-k matrices with at most αk^{β} field multiplications. Here α and β are real numbers with β > 2. Then one can find a quasiinverse of a k-by-k matrix, together with a basis for the nullspace of the matrix, in O(k^{β}) field multiplications.

Theorem 4. Assume that m is a power of 2, that n is a positive multiple of m, and that p is always chosen as m/2 when m ≥ 2 in Algorithm Q. Set γ = 7/(2^{β} − 4). Then Algorithm Q uses at most αγnm^{β−1} field multiplications.

Proof. For m = 1 Algorithm Q performs no multiplications.

For m ≥ 2 the matrix products in Algorithm Q may be broken up into, at worst, n/p + 5 products of p-by-p matrices, using a total of α(n + 5p)p^{β−1} field multiplications. The two recursive calls are each of size at worst p by n, so by induction they take at most 2αγnp^{β−1} field multiplications.

Since p = m/2 ≤ n/2 we have n + 5p + 2γn ≤ (7/2 + 2γ)n = γ2^{β−1}n. □

Under the same assumptions, Algorithm S uses at most α(δ_{1}nm^{β−1}+δ_{2}m^{β})+m^{2} field multiplications, and Algorithm N uses at most α(δ_{3}n^{2}m^{β−2} + (4γ/7)nm^{β−1} + δ_{2}m^{β}) field multiplications. Here δ_{1} = (2 + 2γ)/(2^{β} − 2), δ_{2} = 1/(2^{β} − 1), and δ_{3} = (4 + 2γ)/(2^{β} − 2).

One-dimensional pivots. If p = 1, my algorithms reduce to Gaussian elimination with partial pivoting.

Theorem 5. Assume that p is always chosen as 1 in Algorithm Q. Then Algorithm Q takes at most rb(m, n, r) field multiplications, where b(m, n, r) = mn + r^{2} − 1 − m(r − 1)/2 − n(r + 1)/2.

For r = m this bound factors as m(m − 1)(m + n + 2)/2. For r = m = n the bound is roughly m^{3}; for r = m ≪ n the bound is roughly m^{2}n/2.

Proof. For r = 0 there are no multiplications. For m = 1 and r = 1 there are again no multiplications; and b(1, n, 1) = 0.

For m ≥ 2 we have two cases. If s = 0, Algorithm Q ends up doing nothing except calling itself recursively with an (m − 1)-by-n matrix of rank r; by induction this takes at most rb(m − 1, n, r) ≤ rb(m, n, r) multiplications.

If s = 1, Algorithm Q does various matrix products totalling (m−1)n+2r(r−1) field multiplications. It also calls itself recursively with an (m−1)-by-(n−1) matrix of rank r −1. Total: (r −1)b(m−1, n−1, r −1)+(m−1)n+2r(r −1) = rb(m, n, r) multiplications. □

Similarly, Algorithm S uses at most r(mn−m(r −3)/2−n(r +1)/2+(r^{2} −1)/3) multiplications. For r = m = n this is roughly m^{3}/3. For r = m ≪ n this is roughly m^{2}n/2.

Algorithm N uses at most XXX. For r = m = n this is roughly m^{3}/3. For r = m ≪ n this is roughly mn^{2}.

Practical concerns. XXX Reader beware: I have not yet tested these algorithms; I have merely proven them correct.

How should we choose p in practice? One sensible possibility is p = ⌊m/2⌋. This takes maximum advantage of fast multiplication, for large m. Even without fast multiplication, it is a good idea to break the original matrix into solid chunks rather than peeling off one row at a time, because this uses the cache very effectively.

On the other hand, there is a mild source of inefficiency in Algorithm Q as p grows. The reason is that Algorithm Q always computes all of C in Theorem 2, even though in principle only a portion of C should be necessary. Imagine, for example, that C is a matrix with 100 rows and 10000 columns. Its quasiinverse F is determined by at most 100 of those columns. If we somehow knew which columns to look at, or if we guessed the right columns and verified them against some magically acquired knowledge of C’s rank, we could avoid computing the other columns of C. Note that Strassen’s (failure-prone) method avoids this slight slowdown, essentially by assuming that the first 100 columns of C form a nonsingular matrix. One could perhaps modify Algorithm Q to optimistically check for this case, and to back off if those columns turn out to be singular.

I recommend experimenting with various choices of p between 1 and m/2, for whatever type of matrix is at hand.

Note that there is a dual form of generalized Gaussian elimination, in which V, rather than W, is split up. Perhaps one could gain some speed in practice by combining the original method with its dual.

References

1. J. R. Bunch and J. E. Hopcroft, XXX, Mathematics of Computation 28 (1974), 231–236.

2. Donald Coppersmith and S. Winograd, XXX, Journal of Symbolic Computation 9 (1990), 251–280.

3. Carl Linderholm, Mathematics made difficult, XXX, XXX.

4. Volker Strassen, Gaussian elimination is not optimal, Numerische Mathematik 13 (1969), 354.

X. XXX, XXX, to appear, XXX XXX (XXX), XXX–XXX.

5 Brewster Lane, Bellport, NY 11713