- Email: [email protected]
A generalization of the fast LUP matrix decomposition algorithm and applications
JOURNAL
OF ALGORITHMS
3,45556 (1982)
A Generalization of the Fast LUP Matrix Decomposition Algorithm and Applications OSCAR H. IBARRA* AND SHLOMO MORAN**+ Department of Computer Science, University of Minnesota, Minneapolis, Minnesota AND
ROGERHur* Department of Computer Science, University of Toronto, Toronto, Canada Received May 27, 1981; revised July 25, 1981
We show that any m X n matrix A, over any field, can be written as a product, LSP, of three matrices, where L is a lower triangular matrix with l’s on the main diagonal, S is an m X n matrix which reduces to an upper triangular matrix with nonzero diagonal elements when the zero rows are deleted, and P is an n X n permutation matrix. Moreover, L, S, and P can be found in O(m”-‘n) time, where the complexity of matrix multiplication is O(m“). We use the LSP decomposition to construct fast algorithms for some important matrix problems. In particular, we develop O(m’-‘n) algorithms for the following problems,where A is any m X n matrix: (I) Determine if the system of equations A,? = b (where b is a column vector) has a solution, and if so, find one such solution. (2) Find a generalized inverse, A*, of A (i.e., AA*A = A). (3) Find simultaneously a maximal independent set of rows and a maximal independent set of columns of A.
1. INTRODUCTION
In 1969, Strassen discovered a fast recursive algorithm for multiplying two n X n matrices in O(n2,8’) time [8].’ He also showed that any O(n”) algorithm, where (Y> 2, can be used to obtain O(n”) algorithms for matrix *Research supported in part by NSF Grant MCS78-01736. +Present Address: Department of Computer Science, Technion, Haifa 32OOO,Israel. *Research supported in part by the Natural Sciences and Engineering Research Council of Canada. ‘This bound has since been improved. The best bound known today is due to Coppersmith and Winograd, and is 0( n 2 495.. ).
45 0196-6774/82/010045-12$02OO/O Copyright All rights
8 1982 by Academic Press. Inc. of reproduction in any form resewed
46
IBARRA,
MORAN,
AND
HUI
inversion, computation of determinants, and solution of simultaneous linear equations, provided that each matrix encountered during the recursion process is nonsingular. This last result was later shown to hold for all matrices which are initially nonsingular by Bunch and Hopcroft (see [ 1,4]). Thus, Gaussian elimination is not optimal for solving a system of simultaneous linear equations AZ = b when A is nonsingular. The proof in [l] consisted of showing that any n X n nonsingular matrix A can be decomposed into k = LUP, where L is an n X n lower triangular matrix with l’s on the main diagonal, U is an n X n upper triangular matrix with nonzero diagonal elements, and P is an n X n permutation matrix. Moreover, L, U, and P can be found in O(na) time. (Thus, to solve AZ = b, first solve Ly = 6 for jJ and then solve UPX = jj for 2. Hence, knowing L, U, and P, the solution X can be found in O(n2) time.) When A is singular, an LUP decomposition may not exist, and even if it exists, it may not provide a fast algorithm for solving AZ = b (see Section 2). In this paper, we modify the algorithm of Bunch and Hopcroft [l] to show that any m X n matrix A (m I n) can be decomposed into a product LSP of three matrices, where L is an m X m lower triangular matrix with l’s on the main diagonal, S is an m X n matrix which reduces to an upper triangular matrix with nonzero diagonal elements when the zero rows are deleted, and P is an n X n permutation matrix. Moreover, L, S, and P can be found in O(m a’-‘n) time. We then use this decomposition to develop O(m”- ‘n) algorithms for the following problems, where A is any m X n matrix: (1) Determine if the system AZ = b(where his a column vector) has a solution, and if so find one such solution. This shows that Gaussian elimination is not optimal for solving any system of equations AZ = b: (2) Find a generalized inverse, A*, of A (i.e., AA*A = A). (3) Find simultaneously (in A) a maximal independent set of rows and a maximal independent set of columns. (4) Diagonalize A; i.e., find an m X m nonsingular matrix X and an n X n nonsingular matrix Y such that Jyy=
4 N-1 0
O) 0
where 1, denotes the r X r identity matrix. (5) Find the rank of A. (6) Find the optimal (least square) solution to a (possibly inconsistent) system of equations A5 = 6 over the complex field.
47
FAST LUPMATRIXDECOMPOSITION
2. LSP DECOMPOSITIONOF
MATRICES
Let A be an m X n matrix (m 5 n). An LUP decomposition of A (if it exists) is a triple (L, U, P), where L is a lower triangular m X m matrix with l’s on the main diagonal, U is an upper triangular m X n matrix, P is an n X n permutation matrix, and A = LUP [I]. If A is nonsingular (i.e., rank(d) = m), then an LUP decomposition of A exists, and it can be used to obtain a solution to a system of simultaneous equationsAx = bin O(m*) arithmetic operations. In [l], an algorithm which finds an LUP decomposition of any m X n nonsingular matrix in O(m u-‘n) time is given. This algorithm provides an O(m”-‘n) algorithm to find a solution of AZ = b. However, the algorithm may fail if A is singular, and therefore it cannot be used to decide if a singular system of equations A.Y = b has a solution (and to find one such solution if it exists). Moreover, we observe that if A is singular, then an LUP decomposition of A may not really be useful in constructing a fast algorithm to solve AZ = b. Consider, for example, the system AZ = 6, where
-42nx2n
=[3+],
i=[:‘I. 2n
B is any n X n matrix. A is already upper triangular; hence, A = U = LUP for L = P = Z2”. However, solving fi = his no easier than solving BZ = F, where
Remark. An LUP decomposition of a singular matrix may not always exist. For example, one can easily check that the matrix
does not have an LUP decomposition. Now define an L’U’P decomposition of A, where now L’ is an m X m lower triangular matrix (possibly singular), U’ is an m X n upper triangular matrix with nonzero diagonal elements, and P is a permutation matrix. It can be shown that any m X n matrix A has an L’U’P decomposition, and there is an algorithm to find L’, U’, and P in O(m”‘-‘n) time. However, we
48
IBARRA,
MORAN,
again see that an L’U’P decomposition of equations. For let
A
2nx2n
= [+I0 B
0 0
AND
HUI
may not be useful in solving systems
and
6:
b, II ;
(
b 2n
where B is any n X n matrix. Then A = L’ = L’U’P for u’ = P = Z,,. But then the complexity of solving L’X = b is equivalent to the complexity of solving B.5 = C, where
In this section, we introduce a decomposition called LSP. L and P are as in the LUP decomposition, and S is a matrix which reduces to an upper triangular matrix with nonzero diagonal elements when the zero rows are deleted. When A is nonsingular, LSP reduces to LUP. We shall define this concept formally, and we shall show that an LSP decomposition of any rn X n matrix can be found in O(m”-‘n) time. We shall also show how an LSP decomposition can be used to construct a fast algorithm to decide if Ax = b has a solution, and to find one such solution if it exists. DEFINITION. An m X n matrix S is semi-upper triangular (s.u.t.) if deleting the zero rows from S results in an upper triangular matrix with nonzero diagonal elements. (By convention, the zero matrix [O],.. is s.u.t.) EXAMPLE.
The following matrix is s.u.t. 2
1 0
10
1
s=o
1 0. 2I
DEFINITION. Let S,,,Xn be s.u.t. of rank r > 0. (Clearly, the rank of S is the number of its nonzero rows.) Let S = [S, ] S,], where S, is an m X r matrix. Then an m X m matrix S (-‘) is a left semi-inverse of S if s’-
1)s =
z, .-* L-l--l 00
and
[-1
SC-us, = zr 0
LEMMA 2.1. Let S = [S, ( S,] be s.u.t. of rank r > 0, where S, is an m X r matrix. Then SC-‘) can be computed in O(ma) time.
49
FAST LUPMATRIXDECOMPOSITION
Proof.
The following algorithm
constructs SC-‘):
(1) Permute the rows of S, so that the resulting matrix is of the form u - 7 10 I where U is an r X r upper triangular ments. (2) Compute U-’ and set VmXm
(3) Reverse the permutation resulting matrix is S-l).
=
matrix with nonzero diagonal
u-’ [-cl 0
0 0
done in step ( 1) on the columns of V. The
The correctness of the algorithm is easily verified. Clearly, complexity is O(m + ra + r2) I O(m*).* El EXAMPLE.
ele-
the time
Let 1 0 1
1 0.1 21
Then s, =
II 1 0
3
2
Indeed,
1 1 1
o;-f =o 1,2 -o-&0
and
2The inverse of an r X r triangular matrix can be found in cP + O( r*) steps, where cna steps are required for matrix multiplication [I].
50
IBARRA,
MORAN,
AND
We can now define an LSP decomposition
HUI
formally.
DEFINITION. An LSP decomposition of an m X n matrix A is a triple (L, S, P), where Lx,,, is lower triangular with l’s on the main diagonal, SmXn is s.u.t., and Pnxn is a permutation matrix. We now generalize the algorithm in [ 1, 41 to find an LSP decomposition of an m X n matrix A (m 5 n). Without loss of generality, we assume that m and n are powers of 2. If m = 1 a trivial algorithm of O(n) steps is applied. (If A is a zero vector [0, . . . , 01, then L = [l], S = A, and P = I,.) For a matrix with 2m rows, the following recursive procedure is applied: Let
where A, and A, are m X n matrices. Step 1. Compute L,, S,, P, -the LSP decomposition rank (A,). Compute A; = A, Pr ‘. Then
then let C’ = If s, = K%nXn~
of A,. Let r, =
A, and go to step 5. (In this case, L, = Z, and
P, = Z,.) Otherwise, let
[;I=[% where F and S; are m X r, matrices. Step 2.
Compute S{-‘).
Step 3. Let FLx,,, be the matrix obtained by adding to F m - r, zero columns, i.e., F’ = [F ( 01. Compute F’S{-‘). By Lemma 2.1,
y’s;=I-rrlI. 0
Hence, F’S;-%;
I1
= [ FI 0] 5 0
= F.
51
FASTLUPMATRIXDECOMPOSITION
Step 4. Compute C’ = C - KS,‘-‘)B.
Then
A=[&+$+$. Step 5. Compute L,, S,, P,-the
A=
L, [t 0
0 L2
LSP decomposition
[I 0
P
s,
of C’. If S, = [0],
where L, = I,. 2’
Else
A= L’ H-1 F’S’-” I
0 L
“,; 1 “;‘I[;
1 jP,
= LSP.
2
Time analysis. Let the complexity of n X n matrix multiplication be cn”. (Note that this implies multiplication of D,,,x,,,C,,,xn, where m 1n, takes cm”-‘n steps.) Let I(n) be the complexity of inverting an n X n triangular matrix. (As noted before, I(n) = cn” +‘O(n2).) Then, for m 5 n with m and n powers of 2, we have the following relation for the complexity 7’( m, n) of the algorithm above (we neglect lower-order terms): T( 1, n) I bn for some b (without loss of generality assume b I c). T(2m, n) 5 2T(m, n) + Z(m) + cma + cma-‘n 5 2T(m, n) + 3cm”-‘n. This easily implies (by induction)
that
T(m, K) I 3cm”-In/
(2cr-’ - 2).
From the discussion above, we have our first theorem: THEOREM 2.1. There is an O(m”-’ n) algorithm position of any m X n matrix (m 5 n). COROLLARY
in O(m”-‘n)
2.1. time:
to find an LSP decom-
Let m I n. Each of the following problems can be solved
(1) Given an m X n matrix A, find its rank. (2) Given an m X n matrix A and a column vector 6, determine if the system of equations A% = b has a solution, and if so find one.
52
IBARRA,
MORAN,
AND
Proof. (1) To compute rank(A), rank(A) = number of nonzero rows in ing rank(A) in O(m”-‘n) time is given (2) To solve Ax = 6, decompose A Then solve SPY = jj for X. El
HUI
decompose A into LSP. Then S. (A different method for computin [9].) into LSP. Next, solve LjJ = b for 7.
A variation of the LSP decomposition which is sometimes more convenient to use is given by the following definition. DEFINITION. An LQUP decomposition of an m X n matrix A is a quadruple, (L, Q, U,P>, where Lmxm is lower triangular with l’s on the main diagonal, Q,,, and Pnxn are permutation matrices, and
[1
u= !
0
,
where U, is upper triangular with nonzero diagonal elements. THEOREM 2.2. There is an O(m*-’ n) algorithm to find an LQUP decomposition of any m X n matrix (m I n).
Proof:
Let A = LSP. Let Q, be an m X m permutation
matrix such that
i1
Q,S = ;
,
where U, is upper triangular with nonzero diagonal elements. (By the structure of S, such a Q, exists and can be easily found in 0( mr) steps.) Let Q = Q,-’ and U= Q,S. Then LQUP = LQl’(Q,S)P = LSP. Cl
3. OTHER APPLICATIONS OF LSP/LQUP
DECOMPOSITION
In this section, we use the LSP/LQUP decomposition to provide fast algorithms for the following tasks, where A is an m X n nonzero matrix: (1) Diagonalize A; i.e., find nonsingular
matrices X,,,
and Ynxn such
that Jgy=
L-I--l 00’ 4
O
where r = rank(A). (2) Find a maximal nonsingular square submatrix of A; i.e., find indices 1 5 i,, i,, . . . , i, 5 m, 1 5 j,, j,, . . . , j, 5 n such that the matrix
FAST
LUP
MATRIX
53
DECOMPOSITION
formed by the elements at the intersections of rows i,, . . . , i, and columns j, is nonsingular, and r is maximal. II,..., (3) Find a generalized inverse of A, i.e., a matrix AZ,,, such that AA*A = A [2, 3, 6, 71. Task 1 corresponds to the “classical” matrix diagonalization problem. Task 2 is a generalization of the problem of finding a maximal independent set of vectors (i.e., a “base”) from a given set of vectors. Task 3 has many applications in applied mathematics [2, 3, 6, 71. When A is nonsingular, fast algorithms for these tasks follow rather easily from the LUP decomposition
[Il. First, we consider Task 1. THEOREM
3.1.
There is an 0( ma-’ n) algorithm
to diagonalize an m X n
matrix. Proof:
Let A = LQUP be the decomposition
described in Theorem 2.2,
i.e., U
MXfl
=
‘1 N-1 0
where U, is an r X r nonsingular, and Y = P-‘Y’ 7 where
F =
Q-‘L-‘Ap-1,
0
upper triangular matrix. Let X = Q -‘L-l
Then
XAy=(Jy’=
I--1 00 Ir
’ .
The time to compute X is clearly 0( m”). To compute Y, we first have to compute Y’. Computing U,-*F requires O(r”-‘n) time. Now Y’ has less than n(r + 1) nonzero entries. Thus, P-‘Y’ can be computed in O(nr) additional steps. It follows that diagonalization can be carried out in O(m”-‘n) steps. 0 Note.
The matrix
j.,,= p-1 -u,-‘F I I,-,
54
IBARRA,
MORAN,
AND
HUi
(U, and F defined above) has rank n - r. Moreover, AN = [O]mX(n--r). That is, the columns of N form a basis for the null space of A. Thus, if AZ = b has a solution, then we can find all solutions in the form X0 + NY, where X0 is some solution (obtained, e.g., as in Corollary 2.1) and J is any n - r element vector. We now consider Task 2. We will need the following definition
[5].
DEFINITION. Let A, x n be a given matrix, and 1 5 i,, i,, . . . , i, I m, 1 sj,, j2,. . . , j, I n. Then A[i,, . . . , i, ]j,, . . . , j,] denotes the submatrix of A formed by the elements at the intersections of rows i,, . . . , i, and columns J,,...r j,. If A is of rank r and for some i,, . . . , i,, j,, . . . , j,., A’ = A[i,, . . . , i, li r, . . . , j,.] is nonsingular, then A’ is a maximal nonsingular square submatrix of A.
It is known that for A of rank r, A[ i,, . . . , i, ]j,, . . . , j,] is nonsingular if and only if rows A,,, . . . , Ai, are independent and columns A,,, . . . , Air are independent. The following lemma provides an easy way to find a nonsingular A[i,, . . . , i, ]j,, . . . , j,] from the LSP decomposition of A. LEMMA 3.1. Let A = LSP be of rank r. Let i,, . . . , i, be the nonzero rows P(r) = jr. (P(i) is the image of i under permutation ofS,andP(l)=j,,..., P.)3 Then A[i,, . . . , i, ]j,, . . . , jr] is nonsingular. Proof We need only show that rows A,,, . . . , Ai, are independent, and columns Ajl, . . . , Ail are independent. Clearly, rows S,,, . . . , Si, are independent. One can easily verify that LS[i,, . . . , i, ] 1,2, . . . , n] = L[i,, . . . , i, I J,, . * *, i,] S[i,, . . . , i, I 1,2,. . . , n]. Since L is lower triangular and nonsingular,-so is L[i,, . . . , i, I i,, . . . , i,]. Hence, LS[i,, . . . , i, I 1,2, . . . , n] is nonsingular, which means that rows i,, . . . , i, in LS are linearly independent. Clearly, the multiplication of LS on the right by a permutation matrix P does not change this fact. It follows that A,,, . . . , Ai, are independent. Now columns S’, . . . , S’ are independent, and therefore so are the first r columns in LS. It follows that in A = LSP, columns P(l), . . . , P(r) are independent. n From Lemma 3.1, we have THEOREM 3.2. There is an O(m”-’ n) algorithm to find a maximal nonsingular square submatrix of an m X n matrix.
Finally, we consider Task 3. Recall that A:,,,, is a generalized inverse of A mXn if AA*A = A. A reflexive generalized inverse of Amxn is a matrix A’,,, satisfying AA’A = A and A’AA’ = A’ [3]. In general, a reflexive generalized 3We can consider P as a permutation rather than a permutation matrix.
FAST
LUP
MATRIX
inverse is not unique. However, if A is square and nonsingular, A’ = A-’ and is therefore unique. THEOREM 3.3. There is an O(m”-’ alized inverse of A.
Proof. By Theorem such that
55
DECOMPOSITION
then
n) algorithm to find a reflexive gener-
3.1, we can compute X and Y in time O(m”-‘n)
Then A* is a generalized inverse if and only if
‘x D-1 vw [-cl v vu
A*=y4
for some U, V, W. A’ is a reflexive generalized inverse if and only if
A’=J’I,
’
X
for some U, V (see [2,3]). 0 Let us restrict ourselves now to matrices over the complex field. The Moore-Penrose or pseudoinverse Aix, of A,,, satisfies AA+A = A, A+AA+ = A+, (AA+)”
= AA+,
(A+A)” = A+A, where AH is the conjugate transpose of A. It can be shown that, complex A, At is unique [3].4 A+ has the additional property that X0 = the optimal solution to the (possibly inconsistent) system A? = 6, sense that for all vectors X, II AjTo - ill 5 IIAZ - bll and IIA&, IlAT - &II implies lIZoIl 5 IIXII. THEOREM
for a A+bis in the bll =
3.4. At can be computed in time O(m*-‘n).
PrwoJ Consider A = LQUP. Since the last m - r rows of UP are zeros, A = LQ UP, where @ is the first r columns --of LQ and w is the first r rows of UP. It is easily verified that A+ = vp”( UP UPH) - ‘(aHa) - ‘aH. Cl w1t,,, = PI .x,.) The total time is O(m”-‘n + rQ) = O(m”-‘n). 4At may not exist for matrices over fields with finite characteristic (where AH is interpreted as the transpose of A).
56
IBARRA, MORAN, AND HUI
ACKNOWLEDGMENTS This paper was motivated by a remark made by an anonymous referee that the LLIP decomposition technique in [I] can be used to find the rank of a matrix. (We note that a more direct way to find the rank using the LUP decomposition has been obtained in [9].)
REFERENCES I. A. AHO, J. HOPCROFT,AND J. ULLMAN, “The Design and Analysis of Computer Algorithms,” Addison-Wesley, Reading, Mass., 1974. 2. A. BEN ISRAEL AND T. GREVILLE, “Generalized Inverses, Theory and Applications,” Wiley, New York, 1974. 3. T. BOULLION AND P. ODELL, “Generalized Inverse Matrices,” Wiley-Interscience, New York, 1971. 4. J. BUNCH AND J. HOPCROFT,Triangular factorization and inversion by fast matrix muhiplication, Math. Comp. 28 (1974), 231-236. 5. M. MARCUS AND H. MING, “A Survey of Matrix Theory and Matrix Inequalities,” Allyn & Bacon, Boston, 1974. 6. R. MEYER, “Essential Mathematics for Applied Fields,” Springer-Verlag, New York/Berlin, 1979.
P. PRINGLE AND A. RAYNER,“Generalized Inverse Matrices and Applications to Statistics,” Griffin’s Monograph 28, Hafner, New York, 1971. 8. V. STRASSEN,Gaussian elimination is not optimal, Numer. Math. 13 (1969), 354-356. 9. R. COLE, private communication. 7.
OF ALGORITHMS
3,45556 (1982)
A Generalization of the Fast LUP Matrix Decomposition Algorithm and Applications OSCAR H. IBARRA* AND SHLOMO MORAN**+ Department of Computer Science, University of Minnesota, Minneapolis, Minnesota AND
ROGERHur* Department of Computer Science, University of Toronto, Toronto, Canada Received May 27, 1981; revised July 25, 1981
We show that any m X n matrix A, over any field, can be written as a product, LSP, of three matrices, where L is a lower triangular matrix with l’s on the main diagonal, S is an m X n matrix which reduces to an upper triangular matrix with nonzero diagonal elements when the zero rows are deleted, and P is an n X n permutation matrix. Moreover, L, S, and P can be found in O(m”-‘n) time, where the complexity of matrix multiplication is O(m“). We use the LSP decomposition to construct fast algorithms for some important matrix problems. In particular, we develop O(m’-‘n) algorithms for the following problems,where A is any m X n matrix: (I) Determine if the system of equations A,? = b (where b is a column vector) has a solution, and if so, find one such solution. (2) Find a generalized inverse, A*, of A (i.e., AA*A = A). (3) Find simultaneously a maximal independent set of rows and a maximal independent set of columns of A.
1. INTRODUCTION
In 1969, Strassen discovered a fast recursive algorithm for multiplying two n X n matrices in O(n2,8’) time [8].’ He also showed that any O(n”) algorithm, where (Y> 2, can be used to obtain O(n”) algorithms for matrix *Research supported in part by NSF Grant MCS78-01736. +Present Address: Department of Computer Science, Technion, Haifa 32OOO,Israel. *Research supported in part by the Natural Sciences and Engineering Research Council of Canada. ‘This bound has since been improved. The best bound known today is due to Coppersmith and Winograd, and is 0( n 2 495.. ).
45 0196-6774/82/010045-12$02OO/O Copyright All rights
8 1982 by Academic Press. Inc. of reproduction in any form resewed
46
IBARRA,
MORAN,
AND
HUI
inversion, computation of determinants, and solution of simultaneous linear equations, provided that each matrix encountered during the recursion process is nonsingular. This last result was later shown to hold for all matrices which are initially nonsingular by Bunch and Hopcroft (see [ 1,4]). Thus, Gaussian elimination is not optimal for solving a system of simultaneous linear equations AZ = b when A is nonsingular. The proof in [l] consisted of showing that any n X n nonsingular matrix A can be decomposed into k = LUP, where L is an n X n lower triangular matrix with l’s on the main diagonal, U is an n X n upper triangular matrix with nonzero diagonal elements, and P is an n X n permutation matrix. Moreover, L, U, and P can be found in O(na) time. (Thus, to solve AZ = b, first solve Ly = 6 for jJ and then solve UPX = jj for 2. Hence, knowing L, U, and P, the solution X can be found in O(n2) time.) When A is singular, an LUP decomposition may not exist, and even if it exists, it may not provide a fast algorithm for solving AZ = b (see Section 2). In this paper, we modify the algorithm of Bunch and Hopcroft [l] to show that any m X n matrix A (m I n) can be decomposed into a product LSP of three matrices, where L is an m X m lower triangular matrix with l’s on the main diagonal, S is an m X n matrix which reduces to an upper triangular matrix with nonzero diagonal elements when the zero rows are deleted, and P is an n X n permutation matrix. Moreover, L, S, and P can be found in O(m a’-‘n) time. We then use this decomposition to develop O(m”- ‘n) algorithms for the following problems, where A is any m X n matrix: (1) Determine if the system AZ = b(where his a column vector) has a solution, and if so find one such solution. This shows that Gaussian elimination is not optimal for solving any system of equations AZ = b: (2) Find a generalized inverse, A*, of A (i.e., AA*A = A). (3) Find simultaneously (in A) a maximal independent set of rows and a maximal independent set of columns. (4) Diagonalize A; i.e., find an m X m nonsingular matrix X and an n X n nonsingular matrix Y such that Jyy=
4 N-1 0
O) 0
where 1, denotes the r X r identity matrix. (5) Find the rank of A. (6) Find the optimal (least square) solution to a (possibly inconsistent) system of equations A5 = 6 over the complex field.
47
FAST LUPMATRIXDECOMPOSITION
2. LSP DECOMPOSITIONOF
MATRICES
Let A be an m X n matrix (m 5 n). An LUP decomposition of A (if it exists) is a triple (L, U, P), where L is a lower triangular m X m matrix with l’s on the main diagonal, U is an upper triangular m X n matrix, P is an n X n permutation matrix, and A = LUP [I]. If A is nonsingular (i.e., rank(d) = m), then an LUP decomposition of A exists, and it can be used to obtain a solution to a system of simultaneous equationsAx = bin O(m*) arithmetic operations. In [l], an algorithm which finds an LUP decomposition of any m X n nonsingular matrix in O(m u-‘n) time is given. This algorithm provides an O(m”-‘n) algorithm to find a solution of AZ = b. However, the algorithm may fail if A is singular, and therefore it cannot be used to decide if a singular system of equations A.Y = b has a solution (and to find one such solution if it exists). Moreover, we observe that if A is singular, then an LUP decomposition of A may not really be useful in constructing a fast algorithm to solve AZ = b. Consider, for example, the system AZ = 6, where
-42nx2n
=[3+],
i=[:‘I. 2n
B is any n X n matrix. A is already upper triangular; hence, A = U = LUP for L = P = Z2”. However, solving fi = his no easier than solving BZ = F, where
Remark. An LUP decomposition of a singular matrix may not always exist. For example, one can easily check that the matrix
does not have an LUP decomposition. Now define an L’U’P decomposition of A, where now L’ is an m X m lower triangular matrix (possibly singular), U’ is an m X n upper triangular matrix with nonzero diagonal elements, and P is a permutation matrix. It can be shown that any m X n matrix A has an L’U’P decomposition, and there is an algorithm to find L’, U’, and P in O(m”‘-‘n) time. However, we
48
IBARRA,
MORAN,
again see that an L’U’P decomposition of equations. For let
A
2nx2n
= [+I0 B
0 0
AND
HUI
may not be useful in solving systems
and
6:
b, II ;
(
b 2n
where B is any n X n matrix. Then A = L’ = L’U’P for u’ = P = Z,,. But then the complexity of solving L’X = b is equivalent to the complexity of solving B.5 = C, where
In this section, we introduce a decomposition called LSP. L and P are as in the LUP decomposition, and S is a matrix which reduces to an upper triangular matrix with nonzero diagonal elements when the zero rows are deleted. When A is nonsingular, LSP reduces to LUP. We shall define this concept formally, and we shall show that an LSP decomposition of any rn X n matrix can be found in O(m”-‘n) time. We shall also show how an LSP decomposition can be used to construct a fast algorithm to decide if Ax = b has a solution, and to find one such solution if it exists. DEFINITION. An m X n matrix S is semi-upper triangular (s.u.t.) if deleting the zero rows from S results in an upper triangular matrix with nonzero diagonal elements. (By convention, the zero matrix [O],.. is s.u.t.) EXAMPLE.
The following matrix is s.u.t. 2
1 0
10
1
s=o
1 0. 2I
DEFINITION. Let S,,,Xn be s.u.t. of rank r > 0. (Clearly, the rank of S is the number of its nonzero rows.) Let S = [S, ] S,], where S, is an m X r matrix. Then an m X m matrix S (-‘) is a left semi-inverse of S if s’-
1)s =
z, .-* L-l--l 00
and
[-1
SC-us, = zr 0
LEMMA 2.1. Let S = [S, ( S,] be s.u.t. of rank r > 0, where S, is an m X r matrix. Then SC-‘) can be computed in O(ma) time.
49
FAST LUPMATRIXDECOMPOSITION
Proof.
The following algorithm
constructs SC-‘):
(1) Permute the rows of S, so that the resulting matrix is of the form u - 7 10 I where U is an r X r upper triangular ments. (2) Compute U-’ and set VmXm
(3) Reverse the permutation resulting matrix is S-l).
=
matrix with nonzero diagonal
u-’ [-cl 0
0 0
done in step ( 1) on the columns of V. The
The correctness of the algorithm is easily verified. Clearly, complexity is O(m + ra + r2) I O(m*).* El EXAMPLE.
ele-
the time
Let 1 0 1
1 0.1 21
Then s, =
II 1 0
3
2
Indeed,
1 1 1
o;-f =o 1,2 -o-&0
and
2The inverse of an r X r triangular matrix can be found in cP + O( r*) steps, where cna steps are required for matrix multiplication [I].
50
IBARRA,
MORAN,
AND
We can now define an LSP decomposition
HUI
formally.
DEFINITION. An LSP decomposition of an m X n matrix A is a triple (L, S, P), where Lx,,, is lower triangular with l’s on the main diagonal, SmXn is s.u.t., and Pnxn is a permutation matrix. We now generalize the algorithm in [ 1, 41 to find an LSP decomposition of an m X n matrix A (m 5 n). Without loss of generality, we assume that m and n are powers of 2. If m = 1 a trivial algorithm of O(n) steps is applied. (If A is a zero vector [0, . . . , 01, then L = [l], S = A, and P = I,.) For a matrix with 2m rows, the following recursive procedure is applied: Let
where A, and A, are m X n matrices. Step 1. Compute L,, S,, P, -the LSP decomposition rank (A,). Compute A; = A, Pr ‘. Then
then let C’ = If s, = K%nXn~
of A,. Let r, =
A, and go to step 5. (In this case, L, = Z, and
P, = Z,.) Otherwise, let
[;I=[% where F and S; are m X r, matrices. Step 2.
Compute S{-‘).
Step 3. Let FLx,,, be the matrix obtained by adding to F m - r, zero columns, i.e., F’ = [F ( 01. Compute F’S{-‘). By Lemma 2.1,
y’s;=I-rrlI. 0
Hence, F’S;-%;
I1
= [ FI 0] 5 0
= F.
51
FASTLUPMATRIXDECOMPOSITION
Step 4. Compute C’ = C - KS,‘-‘)B.
Then
A=[&+$+$. Step 5. Compute L,, S,, P,-the
A=
L, [t 0
0 L2
LSP decomposition
[I 0
P
s,
of C’. If S, = [0],
where L, = I,. 2’
Else
A= L’ H-1 F’S’-” I
0 L
“,; 1 “;‘I[;
1 jP,
= LSP.
2
Time analysis. Let the complexity of n X n matrix multiplication be cn”. (Note that this implies multiplication of D,,,x,,,C,,,xn, where m 1n, takes cm”-‘n steps.) Let I(n) be the complexity of inverting an n X n triangular matrix. (As noted before, I(n) = cn” +‘O(n2).) Then, for m 5 n with m and n powers of 2, we have the following relation for the complexity 7’( m, n) of the algorithm above (we neglect lower-order terms): T( 1, n) I bn for some b (without loss of generality assume b I c). T(2m, n) 5 2T(m, n) + Z(m) + cma + cma-‘n 5 2T(m, n) + 3cm”-‘n. This easily implies (by induction)
that
T(m, K) I 3cm”-In/
(2cr-’ - 2).
From the discussion above, we have our first theorem: THEOREM 2.1. There is an O(m”-’ n) algorithm position of any m X n matrix (m 5 n). COROLLARY
in O(m”-‘n)
2.1. time:
to find an LSP decom-
Let m I n. Each of the following problems can be solved
(1) Given an m X n matrix A, find its rank. (2) Given an m X n matrix A and a column vector 6, determine if the system of equations A% = b has a solution, and if so find one.
52
IBARRA,
MORAN,
AND
Proof. (1) To compute rank(A), rank(A) = number of nonzero rows in ing rank(A) in O(m”-‘n) time is given (2) To solve Ax = 6, decompose A Then solve SPY = jj for X. El
HUI
decompose A into LSP. Then S. (A different method for computin [9].) into LSP. Next, solve LjJ = b for 7.
A variation of the LSP decomposition which is sometimes more convenient to use is given by the following definition. DEFINITION. An LQUP decomposition of an m X n matrix A is a quadruple, (L, Q, U,P>, where Lmxm is lower triangular with l’s on the main diagonal, Q,,, and Pnxn are permutation matrices, and
[1
u= !
0
,
where U, is upper triangular with nonzero diagonal elements. THEOREM 2.2. There is an O(m*-’ n) algorithm to find an LQUP decomposition of any m X n matrix (m I n).
Proof:
Let A = LSP. Let Q, be an m X m permutation
matrix such that
i1
Q,S = ;
,
where U, is upper triangular with nonzero diagonal elements. (By the structure of S, such a Q, exists and can be easily found in 0( mr) steps.) Let Q = Q,-’ and U= Q,S. Then LQUP = LQl’(Q,S)P = LSP. Cl
3. OTHER APPLICATIONS OF LSP/LQUP
DECOMPOSITION
In this section, we use the LSP/LQUP decomposition to provide fast algorithms for the following tasks, where A is an m X n nonzero matrix: (1) Diagonalize A; i.e., find nonsingular
matrices X,,,
and Ynxn such
that Jgy=
L-I--l 00’ 4
O
where r = rank(A). (2) Find a maximal nonsingular square submatrix of A; i.e., find indices 1 5 i,, i,, . . . , i, 5 m, 1 5 j,, j,, . . . , j, 5 n such that the matrix
FAST
LUP
MATRIX
53
DECOMPOSITION
formed by the elements at the intersections of rows i,, . . . , i, and columns j, is nonsingular, and r is maximal. II,..., (3) Find a generalized inverse of A, i.e., a matrix AZ,,, such that AA*A = A [2, 3, 6, 71. Task 1 corresponds to the “classical” matrix diagonalization problem. Task 2 is a generalization of the problem of finding a maximal independent set of vectors (i.e., a “base”) from a given set of vectors. Task 3 has many applications in applied mathematics [2, 3, 6, 71. When A is nonsingular, fast algorithms for these tasks follow rather easily from the LUP decomposition
[Il. First, we consider Task 1. THEOREM
3.1.
There is an 0( ma-’ n) algorithm
to diagonalize an m X n
matrix. Proof:
Let A = LQUP be the decomposition
described in Theorem 2.2,
i.e., U
MXfl
=
‘1 N-1 0
where U, is an r X r nonsingular, and Y = P-‘Y’ 7 where
F =
Q-‘L-‘Ap-1,
0
upper triangular matrix. Let X = Q -‘L-l
Then
XAy=(Jy’=
I--1 00 Ir
’ .
The time to compute X is clearly 0( m”). To compute Y, we first have to compute Y’. Computing U,-*F requires O(r”-‘n) time. Now Y’ has less than n(r + 1) nonzero entries. Thus, P-‘Y’ can be computed in O(nr) additional steps. It follows that diagonalization can be carried out in O(m”-‘n) steps. 0 Note.
The matrix
j.,,= p-1 -u,-‘F I I,-,
54
IBARRA,
MORAN,
AND
HUi
(U, and F defined above) has rank n - r. Moreover, AN = [O]mX(n--r). That is, the columns of N form a basis for the null space of A. Thus, if AZ = b has a solution, then we can find all solutions in the form X0 + NY, where X0 is some solution (obtained, e.g., as in Corollary 2.1) and J is any n - r element vector. We now consider Task 2. We will need the following definition
[5].
DEFINITION. Let A, x n be a given matrix, and 1 5 i,, i,, . . . , i, I m, 1 sj,, j2,. . . , j, I n. Then A[i,, . . . , i, ]j,, . . . , j,] denotes the submatrix of A formed by the elements at the intersections of rows i,, . . . , i, and columns J,,...r j,. If A is of rank r and for some i,, . . . , i,, j,, . . . , j,., A’ = A[i,, . . . , i, li r, . . . , j,.] is nonsingular, then A’ is a maximal nonsingular square submatrix of A.
It is known that for A of rank r, A[ i,, . . . , i, ]j,, . . . , j,] is nonsingular if and only if rows A,,, . . . , Ai, are independent and columns A,,, . . . , Air are independent. The following lemma provides an easy way to find a nonsingular A[i,, . . . , i, ]j,, . . . , j,] from the LSP decomposition of A. LEMMA 3.1. Let A = LSP be of rank r. Let i,, . . . , i, be the nonzero rows P(r) = jr. (P(i) is the image of i under permutation ofS,andP(l)=j,,..., P.)3 Then A[i,, . . . , i, ]j,, . . . , jr] is nonsingular. Proof We need only show that rows A,,, . . . , Ai, are independent, and columns Ajl, . . . , Ail are independent. Clearly, rows S,,, . . . , Si, are independent. One can easily verify that LS[i,, . . . , i, ] 1,2, . . . , n] = L[i,, . . . , i, I J,, . * *, i,] S[i,, . . . , i, I 1,2,. . . , n]. Since L is lower triangular and nonsingular,-so is L[i,, . . . , i, I i,, . . . , i,]. Hence, LS[i,, . . . , i, I 1,2, . . . , n] is nonsingular, which means that rows i,, . . . , i, in LS are linearly independent. Clearly, the multiplication of LS on the right by a permutation matrix P does not change this fact. It follows that A,,, . . . , Ai, are independent. Now columns S’, . . . , S’ are independent, and therefore so are the first r columns in LS. It follows that in A = LSP, columns P(l), . . . , P(r) are independent. n From Lemma 3.1, we have THEOREM 3.2. There is an O(m”-’ n) algorithm to find a maximal nonsingular square submatrix of an m X n matrix.
Finally, we consider Task 3. Recall that A:,,,, is a generalized inverse of A mXn if AA*A = A. A reflexive generalized inverse of Amxn is a matrix A’,,, satisfying AA’A = A and A’AA’ = A’ [3]. In general, a reflexive generalized 3We can consider P as a permutation rather than a permutation matrix.
FAST
LUP
MATRIX
inverse is not unique. However, if A is square and nonsingular, A’ = A-’ and is therefore unique. THEOREM 3.3. There is an O(m”-’ alized inverse of A.
Proof. By Theorem such that
55
DECOMPOSITION
then
n) algorithm to find a reflexive gener-
3.1, we can compute X and Y in time O(m”-‘n)
Then A* is a generalized inverse if and only if
‘x D-1 vw [-cl v vu
A*=y4
for some U, V, W. A’ is a reflexive generalized inverse if and only if
A’=J’I,
’
X
for some U, V (see [2,3]). 0 Let us restrict ourselves now to matrices over the complex field. The Moore-Penrose or pseudoinverse Aix, of A,,, satisfies AA+A = A, A+AA+ = A+, (AA+)”
= AA+,
(A+A)” = A+A, where AH is the conjugate transpose of A. It can be shown that, complex A, At is unique [3].4 A+ has the additional property that X0 = the optimal solution to the (possibly inconsistent) system A? = 6, sense that for all vectors X, II AjTo - ill 5 IIAZ - bll and IIA&, IlAT - &II implies lIZoIl 5 IIXII. THEOREM
for a A+bis in the bll =
3.4. At can be computed in time O(m*-‘n).
PrwoJ Consider A = LQUP. Since the last m - r rows of UP are zeros, A = LQ UP, where @ is the first r columns --of LQ and w is the first r rows of UP. It is easily verified that A+ = vp”( UP UPH) - ‘(aHa) - ‘aH. Cl w1t,,, = PI .x,.) The total time is O(m”-‘n + rQ) = O(m”-‘n). 4At may not exist for matrices over fields with finite characteristic (where AH is interpreted as the transpose of A).
56
IBARRA, MORAN, AND HUI
ACKNOWLEDGMENTS This paper was motivated by a remark made by an anonymous referee that the LLIP decomposition technique in [I] can be used to find the rank of a matrix. (We note that a more direct way to find the rank using the LUP decomposition has been obtained in [9].)
REFERENCES I. A. AHO, J. HOPCROFT,AND J. ULLMAN, “The Design and Analysis of Computer Algorithms,” Addison-Wesley, Reading, Mass., 1974. 2. A. BEN ISRAEL AND T. GREVILLE, “Generalized Inverses, Theory and Applications,” Wiley, New York, 1974. 3. T. BOULLION AND P. ODELL, “Generalized Inverse Matrices,” Wiley-Interscience, New York, 1971. 4. J. BUNCH AND J. HOPCROFT,Triangular factorization and inversion by fast matrix muhiplication, Math. Comp. 28 (1974), 231-236. 5. M. MARCUS AND H. MING, “A Survey of Matrix Theory and Matrix Inequalities,” Allyn & Bacon, Boston, 1974. 6. R. MEYER, “Essential Mathematics for Applied Fields,” Springer-Verlag, New York/Berlin, 1979.
P. PRINGLE AND A. RAYNER,“Generalized Inverse Matrices and Applications to Statistics,” Griffin’s Monograph 28, Hafner, New York, 1971. 8. V. STRASSEN,Gaussian elimination is not optimal, Numer. Math. 13 (1969), 354-356. 9. R. COLE, private communication. 7.