T 57.a .1B47 SLCO1 3372 DRC- 12 A r;z't- STRUCTURE ANALYSIS AND PARTIAL REINVERSION OF LARGE JACOBIAN MATRICES IN NONLINEAR PROGRAMMING by Johannes Bisschop and Alexander Meeraus Technical Note No. 12 / Research Project 671-58 August 1979 Development Research Center World Bank 1818 H Street, N.W. Washington, D.C. 20433 The views expressed in this paper are those of the authors and do not necessarily reflect those of the World Bank. STRUCTURE ANALYSIS AND PARTIAL REINVERSION OF LARGE JACOBIAN MATRICES IN NONLINEAR PROGRAMMING by Johannes Bisschop and Alexander Meeraus Development Research Center World Bank 1818 H Street, N.W. Washington, D.C. 20433 August 1979 The views expressed in this paper are those of the authors and do not necessarily reflect those of the World Bank. ABSTRACT This paper concerns itself with efficient inversion and updating procedures for large, structurally equivalent Jacobian matrices as they occur in the generalized reduced gradient method for nonlinear programming problems. It is shown that the Schur complement associated with a diagonal block has a special structure induced by the positions and heights of the spikes inside the block. This induced structure is used to develop strategies for the inversion and updating of Schur complements of varying sizes. The resulting updating procedures do not generate any new nonzero elements, and are more efficient than a complete reinversion of the Jacobian matrix at each iteration of the algorithm. Keywords: Inversion, Updating, Partitioning, Spikes, Bumps, Schur Comple- ment, Jacobian Matrices, Bordered Block Lower Triangular. 1. INTRODUCTION Both the inversion and updating of large sparse matrices have been studied widely as there is an ever increasing need for solving large sparse linear and nonlinear programming problems. Major improvements have occurred over the last 15 years. Among them are the various LU inversion and updating techniques [1], (5], [12] and [14], the concept of supersparseness as a data structure [8], the spike selection algorithms of Hellerman and Rarick [6], [7], and the irreducible block finding algorithm of Tarjan [15]. These concepts have been incorporated into 'most commercial large-scale linear programming codes and some of the experimental nonlinear programming codes [3], [10]. Additional ideas have been generated to suggest further improve- ments. Partitioning methods for the inversion of triangular blocks were suggested by Kevorkian [9] and McBride [11], while a compact updating scheme based on partitioning and augmentation was proposed in [2]. In this paper we would like to extend the use of partitioning techniques for the inversion and updating of sequences of large, structually equivalent Jacobian matrices. Based on our limited experience with solving economic nonlinear optimization models using a sparse implementation of the generalized reduced gradient algorithm [3], [10], we have observed that the zero-nonzero pattern of the Jacobian matrices remains unchanged over an average sequence of at least 5 iterations. In addition, we have observed that most large nonlinear optimization problems contain a sizable proportion of linear terms in the constraint set, which implies that large portions of the Jacobian matrix are invariant from iteration to iteration. These observations naturally lead to the question as to how an implementation of the generalized reduced - 2 - gradient algorithm can take advantage of them. The existing two codes represent opposite extremes. In [10] an LU-factor update is performed for every.basic variable (column) with nonlinear entries in the constraint set. If the number of such "nonlinear columns" is small, the effort spent on updating the inverse Jacobian matrix is minimal. In (3] every Jacobian is reinverted subject to the existing pivot assignment. If the number of nonlinear columns happens to be large, then this updating strategy is certainly more efficient in terms of time and space requirements. A blend- ing of the two extremes is probably a good improvement. We would like to develop a middle ground. On the basis of partitioning information, the location of nonlinear columns, and the nesting of spikes within irreducible diagonal blocks, it is possible to maintain a recursive form of the parti- tioned inverse Jacibian matrix using local updating procedures. In this paper we will develop some new insights into nested parti- tioning schemes, and show how these can be used in large-scale nonlinear programming. The topics vary and may seem unrelated at first sight. They do provide us, however, with a framework for the inversion and updating of large structurally equivalent Jacobian matrices. The next section develops the necessary notation and definitions. In section 3 we show how the nested spike structure (as defined in section 2) within diagonal blocks is carried over to the associated Schur complement resulting in a nested bordered block lower triangular structure. This structure can then be analyzed and rearranged to alter the levels-of nesting. If one wants to establish or re-establish the highest levels of nesting subject to existing triangular assignments within a bump, the procedure described in section 4 can be used to reassign -3 - the spikes. In section 5, we show how structurally equivalent column interchanges in the Jacobian matrix result in rank-1 modifications of the Schur complements. Combining these results in section 6, we examine a set of strategies that can be used for the inversion and updating of the partitioned Jacobian matrix. 2. NOTATION AND DEMINITIONS A sparse matrix can be put into a block lower triangular form such that each diagonal block is irreducible (15]. They are irreducible in the sense that no re-arrangement of rows and columns within each block can split the block into more diagonal blocks. The irreducible block find- ing algorithm of Tarjan is not only efficient, but also guarantees to find the maximum number of diagonal blocks. These blocks vary in size, although our ezperience with real world problems suggests that there are relatively few blocks of size greater than one. Those that are of size greater than one will be referred to as bumps, and tend to be sparse themselves. When- ever bumps are large it is advantageous to rearrange rows and columns such that only few columns contain nonzero elements above the diagonal. Those columns with elements above the diagonal are referred to as spikes. A pictorial example of a set of spikes contained in a single bump is presented in Figure 1. Consider an m x m matrix B representing a typical bump. Number all spikes in B from left to right. Let B contain q spikes indexed by the set I = {1, 2, ..., q}. Each spike i in I corresponds to a column number in B indexed by the set (c, c2, ..., c }. There is also a corresponding row number for the peak of each spike indexed by the set {rl, r2, ..., r } Note that for any spike i , ri < ci by definition. For any two spikes k and Z with ck < c., one and only one of the following conditions will be met: Either (i) rZ < rk ; (ii) rk < r < ck ; or (iii) ck < r. 1* 2 3 4 5 6 . 7 8 9 ........ ........................................ 10 12 13 14 .......... . s ..... ...... 15 16 .......... 17 18 19•• 20" 21 22\ 2 - 27 28\ 29% 30 FIGURE 1: A BUM/IP CONTAIN.'3 2 SPIKES World Bank-15890 - 6 - Pictorially, these correspond to the following situations: C c ck et ck c rkr rF\ rk Lk r. r Case i Case ii Case iii Note that case ii can be made into case i via a simple interchange of spikes k and Z. Definition 2.1: A set of spikes contained in B is said to be properly nested if for any two spikes k and Z with ck < c., one of the fol- lowing inequalities holds: a) r 4 r (Case i) b) ck < r, (Case iii) Spike selection algorithms such as the ones proposed by Hellerman and Rarick [6], [7], always results in a set of properly nested spikes. A stack mechanism is employed for those columns that cannot be assigned as triangular columns. It is used in such a way that the kth spike to enter the stack has always a peak value rk less than or equal to that of any other spike entering the stack subsequently. As the stack is emptied on a last-in first-out basis, only cases i and iii can occur. For ease of presentation we will assume that all spikes within a bump are properly nested. Definition 2.2: Use row and column premutations so that the bump B is partitioned into: S B B B2 B 3 B ,1 - 7 - where B1 to lower triangular and B4 is a square matrix with coluams corres- ponding to spikes only. Then this representation of B is called a triangular partitioning of B. Definition 2.3: Consider a triangular partitioning of B. Then the matrix Q - (B4 - B3B 1B2) is called the Schur complement of B1 in B, or just "Schur complement" for short. The Schur complement plays an important role in the inverse repre- sentation of a partitioned matrix B. Whenever B is nonsingular, the system B B2 'x1 b1 B 3 B4 x 2 b 2 has as its solution x2 . (B34 B3B12 -1 (b2 - B3B 1I 1 x1 = B1 (b1 - B2x2 The pros and cons of this partitioned form of the inverse have been examined by Rose and Bunch [13], and Kevorkian 9]. The following definition will be used in the sequel. Definition 2.4: A matrix S partitioned into s blocks of rows and s blocks of colums is said to be bordered block lower triangular (BBLT) whenever the cells Sij contain zero elements for j > i, j # s. - 8 - 3. NESTED BBLT STUCTURE OF THE SCHUR COMPLEMENTS As we saw in the previous section, the size of the Schur complement corresponding to a triangular partitioning of each bump is equal to the num- ber of spikes. As most bumps observed in real-world problems contain only a few-spikes (usually less than 20), the inversion of each Schur complement is a straightforward matter. There are occasional problems that produce bumps with many spikes, so that the corresponding Schur complements are large matrices. One such example is the Indus Basin Model [4]. This model is a linear program- ming problem with 7314 equations and 9616 variables. From observations, we found that a typical basis of this model contains 1 bump of siz& 2700 with roughly 500 spikes. If the corresponding Schur complement happens to be a "full" matrix (a matrix with nonzero elements only), its inversion would be prohibitive in terms of time and storage. Based on our analysis, it is not a full matrix, and has a special zero-nonzero structure of its own. The follow- ing theorem is fundamental in the identification of this structure, and plays an important role in the determination of inversion strategies for large Schur complements. In essence it states that the structure induced by the positions andheightsof the spikes inside a bump, is carried over to the Schur complement. Theorem 3.1: Consider an m x m bump B and a set I of properly nested spikes i, together with their original colum position c. and their peak r., 3. i = 1, 2, ..., q . Consider also the Schur complement Q = (B4 - B3 B 1-1B 2) associated with the triangular partitioning of B. Then for any two spikes k and Z with ck < r., the (k, 9)th entry of Q is zero. th Proof: The (k, c) entry of the matrix Bis the same as the (ck, c2.)t entry of B prior to the triangular partitioning of B. Since ck < r - 9 - (Case iii in the previous section), the (ck, c)th entry of B is zero by assumption. We only need to show, therefore, that the (k, Z)th entry of B3B 1 is also zero. Premultiplying B1-1 with the kth row of B3 and post- multiplytag the result with the Zth column of B2 will produce the (k, Z)th entry of B 3 1 2. Note that the kth row of B3 corresponds to the original c kth row of B with all its spike-related elements removed. Therefore, by construction, the kth row of B3 has all zero's in its last [n-ck-(q-k)] entries. Similarly, the Zth column of B2 corresponds to the original c th column of B with all its spike-related elements removed. Therefore, by construction, and using the assumption that Ck < r, the Zth column of B2 has all zero's in at least the first fck-k] of its entries. Note that the inner product of the kth row of B3 and the Zth colum of B2 is zero, as the total number of consecutive zero elements starting at opposite ends in each of these two (n-q) dimensional vectors is at least (n-q). At this point we merely need to note that, since B1-1 is lower triangular, all consecutive zero elements at the beginning of the Zth column of B2 will be retained after this column has been transformed by B1-. This concludes the theorem. The zero-nonzero element structure induced by Theorem 3.1 is displayed in Figure 2 for the bump of Figure 1. It illustrates how the nested structure of the spikes is carried over in its entirety-to-the Schur complement of the bump. The result is a nested bordered block lower triangular (BBLT) structure as is emphasized by the dotted lines. A "level one" BELT structure always refers to the entire Q-matrix, and represents the partitioning into diagonal blocks and a single border. In Figure 2, the level one diagonal blocks start in columns 1, 5 and 10 respectively, while the border is comprised of columns 29 and 30. A level two BBLT structure always refers 1 2 3 4 5 6 7 8 9 10 1 12 13 14 15 16 17 18 19 20 21 223 24 25 26 27 2829 30 2 3 4 6 % tt 7 8 9 10 12 13 14 [ 15 16 i 1 2 2 W9 17 18 19 .y - 20 21 W 22 23 L 24 25 - 27. 28 29ý .30 warW 8an,-is~es ZERO ELEMENT. FIGURE 2: THE STRUCTURE OF THE MATRIX Q FOR THE BUMP IN FIGURE 1 - 10 - to a level one diagonal, block which is also partitioned into diagonal blocks and a border. Colmns- 17.through 23 in Figure 2 are an example of a level four BBLT structure contained in a level three diagonal block comprised of colums 15 through 27. The arrangement of the spikes in Figure 2 is not the only represen- tation of properly nested spikes. One can substitute out higher levels of nesting in order to obtain more diagonal blocks at a lower level. As a result, the border at this lower level will grow in size. For example, substitute out all levels greater than one in Figure 2. This results in a properly nested level one BBLT structure having 10 diagonal blocks and a border of 6 colunis. The border will contain (from left to right) columns 23, 27, 28, 9, 29, and 30. As we noted already in section 2, the Hellerman-Rarick routines empty the spike stack on a last-in first-out basis, thereby creating the maximum number of nested levels of the BBLT structure. Once we know this arrangement, we have all the flexibility we need for finding alternative arrangements with fewer levels of nesting. Since we will use the nested BBLT structure as a skeleton for the inverse representation of the Schur complement, we have developed in the next section a simple procedure which can recover the maximum number of nested levels whenever this need arises. - 12 - 4. SPIE INSERTIONS WITH F=XED TRIANGULAR ASSIGNMENTS This section developes an algorithm to identify the maximum number of nested spike structures subject to an existing triangular assignment for part of the rows and columns in the Jacobian matrix. There are instances In which such an algorithm is of interest. Whenever tempo- rary spike interchanges have been made for numerical reasons, one can restore the original nesting assuming no structural changes took place. Even if a few structural changes did take place inside the Jacobian matrix, the-method-has merits. Rather than trying to discover the new underlying structure using the relatively time consuming partitioned preassigned pivot procedure of Hellerman and Rarick, one can obtain a good approximation of this structure by leaving all triangular assignments fixed, and inserting the spike columns (rows) at the appropriate places in the triangle. While describing the procedure, we will use the 5 x 5 bump contained in the 12 x 12 sample matrix of Hellerman and Rarick as an example [6]. The column numbers correspond to the numbering of the original matrix. 2 4 8 5 7 2 X X X 3 X X X X B1 B2 8 - X X 6 X X I X 4 X. I X Here we assume the triangular assignment (2, 2) and (3, 4). All remaining colins and rows are as yet unassigned. - 13 - Step 1. Order the unassigned rows according to the number of consecutive zero elements found in the last entries of B3 (in decreasing order). Th-s number determines the position in the triangular part after which the correspondtng row can be inserted. For the above example, we get the following picture. 2 4 8 5 7 2 X X X 3 X X 4 X x 8 X X X 6 X X Row 4 has two consecutive zeros, and can be inserted anywhere after posi- tion 0. Row 8 has one consecutive zero, and can be inserted anywhere after the first triangular row. Row 6 with no zero element at the end can only be placed at the conclusion of the triangular section. Step 2. Consider the unassigned rows in the order determined by step 1, and insert them in their earliest positions. This gives us the set {c, c2, ..., c } described in section 2. The next picture illustrates step 2 for the example. 2 4 8 5 7 4 - X X 2 x X x 8 x x x 3 x x x x 6 x x - 14 - In this case {cl, c2, c3 = {1, 3, 5}. Step 3. Determine the peak values r . Consider the positions c in order, and insert spike k with the largest peak value rk< c. in position ci. This will generate the set {r.} described in section 2. 1 Note that when rk=c, we have a triangular assignment, and the triangular portion of the matrix has grown by one row and one column. Continuing our 5 x 5 example, we are faced with a tie when row 4 is inserted at the beginning of the triangular portion of the matrix. Depending on how the tie was broken, one of the following configurations will result after we complete step 3. 7 2 8 4 5 5 2 8 4 7 4 XJ X X 2 X X X X X 8 X X X X X X 3 X X X X X X X X 6 X X Note that the triangular portion of the matrix has grown by one row (row 4) and one column (either column 7 or column 5). Even though the above three steps guarantee a nested spike structure,- a slight modification of steps 2 and 3 may increase the number of new triangular assignments. Consider the following illustration: - 1.5 - etc. B1 B2 X 0 00 XXXX X 0 0 0 XX X 0 0 0 XX B 3B etc. If the rows in D are inserted in their stated order, the triangular portion B will only grow by 1 row and 1 column. One can verify that if the first row is not inserted first, the triangular portion will grow by 2 rows and 2 coluin. This observation leads to the following modification of steps 2 and 3 Steps 2 and 3 extended. If there are several rows to be inserted in position c., examine the submatrix in B that'is the intersection of these rows and the columns k with current peak value rk > ci. Use the preassigned pivot procedure of Rellerman and Rarick on this rectangular subblock to obtain a triangular assignment involving, say, p rows and p columns. Insert these p rows and columns in the order induced by the triangular assignment starting at posi- tion c . The remaining rows, if any, can then be inserted in any order with spike columns selected according to the usual criterion. As a result of this modification, the number of new triangular assignments will tend to increase. This gain should be weighed against the effort required to identify, extract and analyze the submatrices of B . Once we have identified the sets {c1, c2, ..., c I and {rl, r2, ..., r } for the q spikes, we have all the information we need to analyze the nested BBLT structure, and to substitute out any levels of nesting, if this is desired. - 16 - 5. RANK-ONE MODIFICATIONS WITHOUT STRUCTURAL CHANGES. Whenever two columns have nonzero entries in exactly the same set of rows, they are said to be structurally equivalent, Any time a Jacob±an matrix is evaluated at a new point, certain entries La the non- linear columns will change their numerical value. One can view these changes in the matrix as a sequence of exchanges involving structually equivalent columns. In this section we would like to study the effect of these exchanges an the Schur complement. There are three cases to consider. The simplest case is that of a triangular column which is not part of any bump. If such a column is interchanged for a structurally equivalent column, no additional work is needed to update the Jacobian inverse representation (we defer any numerical problems until the next section). The other two cases are that either a spike or a triangular column within a bump are replaced with a structurally equivalent column. As we shall show, every column interchange within a bump always results in a rank-1 modification of the Schur complement Q. To verify that a spike interchange within a partitioned bump results in a rank-1 modification of Q is straightforward. Let qj be the jth column of Q, where j refers to the position of the spike. Also, let q (b -3T 2), where b and b denote the modified columns in B and j 4 31 2 2 4 2 B respectively. In addition, let v - q q,, and w - e , where e is a unit sector with the one in position j. Then the new Schur complement ^ t Q can be written as Q - Q + vw , which is a rank-1 modification of Q. To show that a triangular column interchange within a bump also results in a rank-1 modification of Q is less straightforward. Consider the . 17 - following notation. Let b denote the newly inserted column in Bl, while b1 denotes the original vector removed from B1. Similarly, let b denote the newly inserted column in B3 and let b3 be the original vector removed from B Then the new Schur complement Q can be written as Q = B4 - [B3 + T T -1 (b3 - b3) e ] B + (b1 - bl)e ] B2 , where e. is a unit vector of appro- priate length, and j refers to the position of the modified column in B1. Using the Sherman-Morrison formula (see e.g. [13]), we may write [B1 + (b1 T -1 -1 -1^ T -1 T -1 b)e ] = B 1 1 (b - bl)e 1 , where S = 1 / (1 + e B 1 b T T -1^-^ bl)) = - (b lj lb j. Let w = 1 B 2 , and let v = 6 (b3 - B3B 1bl Then it is straightforward to verify that the new Schur complement can be a T written as Q = Q + vw , which is a rank-1 modification of Q. As we know from Theorem 3.1, the Schur complement Q has a nested BBLT structure induced by the height and position of each spike. As the spike structure is not affected at all by structurally equivalent column interchanges within a bump, it must be that the modified matrix Q has the same nested structure as Q. This observation implies that the above rank-1 modifications must have a special structure themselves. In the case of a spike interchange within a partitioned bump, one can note that the height of q is the same as that of q by construction. The case of a triangular column interchange within a bump is not so clear cut. One can show for this case that the above modification vw is such that the kth element of w is zero whenever spike k, prior to the triangular partitioning of the bump B, had a peak value rk > j (j refers to the original position of the triangular column), and that - 18 - the kth element of the vector v is zero whenever spike k, prior to the triangular partitioning of B, had a column value ck < j . As an illustration consider the aterchange of a triangular column somewhere between spikes 15 and 16 tn Figure 1. This produces the following zero- TT nonzero pattern for the vectors v and w in the modification vwT. v T: (0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,X,X,X,X,X,X,X,X,X,X,X,X,X,X,X). T wT : (X,X,X,X,X,X,X,X,X,X,X,X,X,X,X,0,0,0,0,0,0,0,0,0,0,0,X,X,X,X). It is of interest to point out that the above rank-1 modifications can be applied to a block recursive structure. Let us return again to Figure 1. Assume that we have a 5 x 5 Schur complement for the block with spikes 5 through 9 as its border, and a 2 x 2 Schur complement for the entire bump with spikes 29 and 30 as its border. If we want to exchange a triangular column from within the inner block with a new but structurally equivalent column, we have to modify both the 5 x 5 and the 2 x 2 Schur complements. This requires two levels of partitioning, implying the -existence of a lower triangular B1 matrix for the inner block, say B11, and a block lower triangular B 1 matrix for the outer block, say BB . An inspection of the modification formulas show that the 2 x 2 complement can be updated using the existing inverse representation for BB , and that the 5 x 5 complement can be modified using BI . Only after the modifications have been performed does the actual interchange of the structurally equivalent columns need to take place. We will return to the notion of nested Schur complements in the next section. - 19 - 6. UZVERSION AND UPDATING STRATEGIES FOR STUCTUALLY EQUIVALENT JACOBIAN MATRICES, In thi1s section we will distinguish between bumps with a few spikes, say less than 20, bumps with a medium number of spikes, say anywhere between 20 and 80, and bumps with many spikes, say more than 80. The choice of these numbers is somewhat arbitrary, but was used to provide some intuitive notion about small, medium and large. Based on our experience with a variety of economic planning models, bumps usually have fewer than 20 spikes. A triangular partitioning of such bumps will produce a relatively small Schur complement Q which can be kept in factored form (either orthogonal or Lower-Upper). As the overall Jacobian matrix is re-evaluated from iteration to iteration, some elements in the bump may change. If there are few nonlinear columns with changing entries inside the bump, one can update the factors of Q using a sequence of rank-1 modifica- tions. If there are many such nonlinear columns, it is better to reinvert the Schur complement. In either case, there is no growth in terms of newly generated nonzero elements. Based on operation counts, the break-even point is when the number of nonlinear columns with changing entries within the bump is roughly one-third the number of spikes. The resulting inversion and updating scheme for these bumps is efficient in terms of storage. As all spikes are grouped together, the partitioning information required for acces- sing the partitioned form of the inverse is kept to a minimum. As the itera- tions progress, it is possible that the Jacobian matrix temporarily drops rank or becomes near-singular. In such a case, one usually inserts one or more logicals to avoid numerical problems. This requires the removal of rows - 20 - and columns. As most data structures for sparse %ystems provide efficient access-either via raws or via columns but seldomly both, the augmentation scheme proposed in t2] can be used to add the logicals. The existing par- t1t1oning of columns and rows remains essentially unchanged as singular rows and columns stay frozen tn the inverse representation. Next we would like to comment on bumps with a medium number of spikes. As the inversion (factorization) of each Schur complement is of the order q3 (q refers to the size of the complement), it becomes advanta- geous to consider the nested BBLT structure of Q. Rather than factorizing the entire Q-matrix, one can limit the inversion to diagonal blocks and associated Schur complements. If the induced BBLT structure produces singu- lar blocks, it is necessary to make minor changes in the block structure. These changes increase the size of the border, and are usually determined during the inversion process. The entire inverse representation of Q can be viewed as a block product form of the inverse. Although there is an improvement in computational time for the generation of this form of the inverse, there is a cost associated with maintaining the relevant partition- ing information for the matrix Q. The updating of the inverse representation under structurally equivalent column interchanges is once again a function of the number of nonlinear columns in the bump. The computation of the break-even point is more complicated in this case as it depends on the sizes of the blocks to be inverted. Observing that the inverse representation is built much faster when only portions of Q are factored, one would expect that on the average the break-even point is when the number of nonlinear columns is much less than one-third of the number of spikes. - 21 - The remaining portion of this section is devoted to bumps with many spikea. Although these bumps are a rare occurrence in reality, they do exist as we observed from the Indus Basin model. Besides containing many sptke*, these bumps are usually large in size, As a result any exist- ing product form or LU inversion scheme will generate a large number of new nonzero elements in the spike columns. Unfortunately this is also true for the above implementations of the partitioned form of the inverse, as the associated Schur complement is large. With the exception of the induced zero elements (see Theorem 3.1), this matrix tends to be full. That is why we propose a different implementation of the partitioned inverse for bumps with many spikes, Rather than storing all the nonzeros of the large Schur complement, it is possible to regenerate them as needed by breaking up the bump. Using the nested BBLT structure as a guide, one can identify bumps within bumps, and compute a separate partitioned inverse representation for each of them. Some minor rearrangements of inner bumps may be required if they turn out to be singular or near-singular, resulting in some growth of outer borders. Referring once again to Figure 1, one can visualize 3 inner bumps and an associated 2 x 2 Schur complement for the outer bump. The third inner bump can in turn be viewed as consisting of 3 inner bumps and a 1 x 1 Schur complement. By inverting only inner bumps we generate relatively few nonzero elements. There is a definite cost, however, associated with break- ing up large bumps. First of all, partitioning information involving several levels must now be kept for all the columns in the bump. Secondly, there is an increased computational cost. If one examines the formulas for the parti- -1 tioned form of the inverse in section 2, one observes that the B matrix 1 needs to be accessed twice. If a bump is nested within an outer bump, the - 22 - B1-1 matrix for the inner bump will be. accessed four times. If the nesting 1_ goes k levels deep, the most iner B 1matrix wll be accessed 2 times. One can remedy thts frequent access by% storing at same level(s) the trans- -1 formed values B 1 B This, of course, generates new nonzero elements. Saving these elements turns out to be equivalent to storing portions of the large Schur complement obtained when all spikes within the bump are grouped together. Decisions as to which transformed elements should be saved, is mostly determined by available storage and the level of nestings employed. Updating the recursive bump structure is again a function of the number of nonlinear columns affecting the bump. The inner bumps are usually small which increases the likelihood that they will be unaffected by the nonlinear columns. Those Schur complements that are affected by the nonlinear columns can either be modified or completely reinverted. The break-even point is the same as the one discussed at the beginning of this section for small Q-matrices. Although our treatment of the above partitioning schemes can at best be viewed as a skeleton for code development, the ideas presented in this paper should be of interest to anyone working on sparse matrix inversion and updating techniques. -23 - REFERENCES [1] Bartels, R.H. and G.H. Golub, "The simplex method of linear programming using LU decomposition". Comunications of the Association for Computing Machinery 12 (1969) pp. 266-268. [2] Bisschop, J. and A. Meeraus, "Matrix augmentation and partitioning in the updating of the basis inverse", Mathematical Programming 13 (1977) pp. 241-254. [31 Drud, A., "An optimization code for nonlinear econometric models based on sparse matrix techniques and reduced gradients", Annals of Economic and Social Measurements 6 (1978) pp. 563-580. [4] Dul6y, J.H., and G.T. O'Mara;~ S.S. Ting and A. Brooke, "Programming and designing investment: the Indus Basin", RPO 671-45, Development Research Center, World Bank, June 1979. [5] Forrest, J.J.H. and J.A. Tomlin, "Updated triangular factors of the basis to maintain sparsity in the product form simplex method", Mathematical Programming 2 (1972) pp. 263-278. [6] Hellerman, E., and D. Rarick, "Reinversion with the preassigned pivot procedure", Mathematical Programming 2 (1971) pp. 195-216. [7] Hellerman, E., and D. Rarick, "The partitioned preassigned pivot pro- cedure"., in D.J. Rose and R.A. Willoughby, eds., Sparse matrices and their applications (Plenum Press, 1972). [8] Kalan, J .E., "Aspects of large-scale in-core linear programming", Proceedings of the 1971 Annual Conference of the Association for Gomputing Machinery, (Chicago, Illinois, 1971). (9] Kevorkian, A.K., "A decompositional algorithm for the solution of large systems of linear algebraic equations", Proceedings of the 1975 IEEE International Symposium on Circuits and Systems, (Boston, Massachusetts, 1975). [10] Lasdon, L.S., and A.D. Warren, "Generalized reduced gradient software for linearly and nonlinearly constrained problems", in H. Greenberg, editor, Design and Implementation of Optimization Software (Sijthoff and Noordhoff Publishers, 1978). [11] McBride, R.D., "A bump triangular dynamic factorization algorithm for the simplex method", Working Paper No. 19, University of Southern California, December 1977. - 24 - [12] Reid, J.K., "A sparsity-exploiting variant of the Bartels-Golub Decomposition for linear programming bases", Report CSS, Atomic Energy Research Establishment (Harwell, Didcot, Oxford- shire, 1975). [13] Rose, D.J., and J.R. Bunch, "The role of partitioning in the numerical solution of sparse systems", in D.J. Rose and R.A. Willoughby, eds., Sparse matrices and their applications (Plenum Press, 1972). [14] Saunders, M.A., "The complexity of LU updating in the simplex method", in R.S. Anderssen and R.P. Brent, eds., The complexity of computational problem solving (University Press, Queensland, 1972). [15] Tarjan, R.E., "Depth first search and linear graph algorithms", SIAM Journal of Computing 1 (1972) pp. 146-160.