UNIVERSITY OF, 
 
 ILLINOIS LIBRARY 
 
 A" URBANA-CHAMPAIGN 
 
 ENGINEERING 
 
NOTICE: Return or renew all Libra 
 each Lost Book is $50.00. 
 
 i Minimum Fee for 
 
 an Materials! Jhe Ulnlmu 
 
 Jul 6 6 1988 
 
 The person charging this material is responsible for 
 its return to the library from which it was withdrawn 
 on or before the Latest Date stamped below. 
 
 Theft, jm^tigrv, aria uriderlining of books are reasons for discipli- 
 ne in 4ismtesatfrorrv the University. 
 Center, '3^H»4pO/ 
 
 nary a 
 To ren 
 
 UNIVE 
 
 \NA-CHAMPAIGN 
 
Digitized by the Internet Archive 
 
 in 2012 with funding from 
 
 University of Illinois Urbana-Champaign 
 
 http://archive.org/details/efficientalgoritOOgilm 
 
Center for Advanced ComDutatior 
 
 UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN 
 
 URBANA. ILLINOIS 61801 /* 
 
 224-2 
 
 CAC Document No. 235 
 
 AN EFFICIENT ALGORITHM FOR SOLVING THE 
 LINEAR INPUT OUTPUT EQUATION WITH 
 AN EXTENSION TO THE NONLINEAR 
 INPUT OUTPUT MODEL 
 
 by 
 Doug Gilmore 
 
 August, 1977 
 
 *»E UBHARY Of IHE 
 
CAC Document No. 235 
 
 An Efficient Algorithm for Solving 
 the Linear Input Output Equation with 
 an Extension to the Nonlinear Input Output Model 
 
 by 
 Doug Gilmore 
 
 Center for Advanced Computation 
 University of Illinois at Urb ana- Champaign 
 Urbana, Illinois 6l801 
 
 August, 1977 
 
 This work was supported in part by the National Science Foundation under 
 Contract US NSF C-10*+5j and was submitted in partial fulfillment of the 
 requirements for the degree of Master of Science in Electrical Engineering 
 in the Graduate College of the University of Illinois at Urbana-Champaign, 
 
 1977. 
 
Ill 
 
 ACKNOWLEDGEMENTS 
 
 I would like to express my appreciation first to Profes- 
 sor Hugh Folk for his approval of the work done on the subject of 
 this thesis. At the beginning many thought that performing such 
 a large computation on a small computer installation was absurd, 
 and I appreciate that Professor Folk gave me the opportunity to 
 demonstrate that the contrary was quite true. I would also like 
 to thank Professor Achmed Sameh for his assistance during the 
 development of the algorithms, and for his comments during the 
 drafting of the thesis. I also thank my thesis advisor, Profes- 
 sor William Perkins for his comments on the drafts of the thesis. 
 His comments have greatly improved the readability of the 
 mathematical text. I also give special thanks to Steve Holmgren. 
 He was unceasingly assistive during the coding of the algorithms. 
 Though I value the assistance of all the above, I alone accept 
 responsibility for any errors. 
 
IV 
 
 TABLE OF CONTENTS 
 
 Introduction 1 
 
 Chapter 1 - Factorization Modification 7 
 
 1.1 - Matrix Decompostion by Elementary 
 
 Transformations and Simplifications thereof 
 
 on 10 matrices 7 
 
 1.2 - Multiple Row and Column Modification 15 
 
 Chapter 2 - A Sensitivity Transformation for Studying the 
 the Perturbations of Solutions Due to 
 Perturbations of the System 21 
 
 2.1 - Perturbations of a Sinqle Element 21 
 
 2.2 - Column Perturbations 28 
 
 2.3 - Row Perturbations 32 
 
 Chapter 3 - A Nonlinear Input Output Model 37 
 
 Discussion and Conclusions 50 
 
 References 53 
 
INTRODUCTION 
 
 Due to recent economic developments such as the energy 
 crisis and persistent cost push inflation there has been consid- 
 erable interest or the structure of national economies, and how 
 crisis effects economic activities. One method that is used to 
 study economic systems is Input Output (10) analysis which is 
 particularly powerful in the study of interindustry activity. 
 The purpose of this paper is to discuss the solution of the 10 
 equation using several methods, the goal being to develop a 
 method to solve a nonlinear Input Output model which will partly 
 alleviate some of the restrictions of linear Input Output 
 analysis. 
 
 Input Output analysis is an econometric method which at- 
 tempts to explain all industrial activity by a simple cause ef- 
 fect relationship. The first critical assumption of 10 analysis 
 is that all goods in a product group are manufactured in an 
 identical manner. From this point on when the term good is men- 
 tioned, it is to be taken as some representative good in one of 
 the product groups. The amount of a good that society needs to 
 produce is the amount to be supplied for final demand plus the 
 
 By final demand we mean several things: goods 
 consumed by households and government, goods sold for 
 export, and also goods used for investment. Most logically 
 investment would be treated as an input to production, but 
 investment can be a very nonlinear function of output. Thus 
 in general economists find it easier to determine investment 
 demands exogeneously . It is possible that the 
 representation ot investment as a nonlinear function of 
 output could make the endogenous determination of investment 
 demands more realistic. 
 
amount used in the production of other goods. The second criti- 
 cal assumption of 10 analysis is the following: when a good is 
 used in the production of another good, the amount used in this 
 activity is always in a fixed proportion to the total production 
 of the other good. These proportionality constants are termed 
 the technical coefficients, and since in general it is possible 
 that all goods can be used directly in the production of all oth- 
 er goods, if we have an n good economy then there are n* techni- 
 cal coefficients. Vvhen the technical coefficients are arranged 
 in a nxn matrix, this matrix is called the technical coefficients 
 matrix or A matrix, though this should not be confused with the A 
 matrix found in control systems literature. If our unit of value 
 is dollars then the element a. . represent the dollar value of the 
 i good required in the direct production of one dollar of the 
 j good. Thus the a^'s are positive fractions. The sum of the 
 
 t" h 
 
 j column of the A matrix represents the fraction of the total 
 cost of producing the good j that is embodied in the goods used 
 in direct production of the j good. Then 
 
 v, = 1 - I a-, (0.1) 
 
 J i=i *J 
 
 represents the unit cost of production not embodied in the goods 
 used in the direct production of the j good. V. is termed the 
 value added for good j and it consists of labor costs, interest, 
 (on the capital used in the direct production of the j tn good,) 
 direct business taxes, and profits. 
 
If y^ is the amount of the i good sold to final demand, 
 then the equation which represents the total production of the 
 i good, that is, x., is 
 
 *i s y t «■ 2 a..x, . (0.2) 
 
 The summation term in (0.2) is the amount of good i needed in 
 the direct production of all goods. Proceeding in the same way 
 for the other n-1 goods, the total demand for all goods is the 
 solution of the matrix equation 
 
 Ax + y = x or (I - A)x = y . (0.3) 
 
 We reiterate the basic assumptions of 10 analysis; 1) that all 
 goods in the same product class are assumed to be made in the 
 same way, 2) that the amount of an input good used in the direct 
 production of another good is always in a fixed proportion. 
 
 The first question that needs to be answered is whether 
 solutions to (0.3) do indeed exist. But since we also desire 
 
 that solutions must correspond to actual economic behavior, the 
 
 2 
 solution vectors should be positive if the final demand vector 
 
 is positive. Bellman [1, pp.296, Theorem 6] has shown that if the 
 
 column sums of A are less that 1, then there is a unique positive 
 
 o 
 
 A vector is said to be positive if all of its 
 
 elements are positive. 
 
solution vector for each positive right hand side. Furthermore, 
 the eigenvalues of A lie inside the unit circle, and 
 
 (I - A)" 1 = I + 2 A 1 , (0.4) 
 
 i = l 
 
 with lim A m = 0, e.g. see Isaacson and Keller [15, pp.15, Theorem 
 
 nHro 
 5]. If k terms are used to approximate (I - A)" 1 then (k-l)n 
 
 multiplications must be performed. 
 
 In the past, eg. Chenary and Clark[17], (I-A) was crude- 
 ly approximated by this method. The objective of this paper is 
 to discuss a method by which the 10 equation can be solved with 
 far more efficiency and accuracy. By a factorization method 
 which will be discussed in Chapter 1, only J:n^ multiplications 
 are required to completely factorize the matrix into a product of 
 
 triangular matrices, which can then be solved for arbitrary right 
 
 2 
 hand sides in n multiplications. The first section of Chapter 1 
 
 will also analyze the numerical stability properties of factori- 
 zation of Input Output matrices by elementary transformations, 
 (sometimes denoted as LU decomposition,) and will reexamine the 
 property of the A matrix so that (0.3) will admit only positive 
 solutions for arbitrary positive final demand vectors. The 
 second section in this chapter proposes a simple method which 
 
 3 
 
 It is standard practice in numerical analysis to only 
 
 count the number of multiplications or divisions in a 
 computation since they are more time consuming than addition 
 or subtraction operations. 
 
will allow modifications of particular rows and columns of the A 
 matrix that will not require the complete ref actor ization of the 
 matrix. This method has a very useful property that allows the 
 successive updating of the 10 matrix without any algorithmic com- 
 plexity. Procedures of this type are generally referred to as 
 factorization modification, and several papers have been written 
 on this subject, including Gill and Murray[13], Golub et al.[12], 
 Sameh and Bezdek[2], and Noh and Sameh[3]. But with the excep- 
 tion of Gill and Murray[13], these only deal with factorization 
 by orthogonal transformations, which will be more numerically 
 stable for general matrices, but require more storage and compu- 
 tation than factorization by elementary transformations thus mak- 
 ing them comparably more expensive. By examining the structure 
 of 10 matrices, we shall see that the use of elementary transfor- 
 mations in the factorization of 10 matrices is quite appropriate. 
 We do not mean to distract the reader away from orthogonal 
 transformation factorization methods. As Bierman[5] points out, 
 (this article is an excellent survey of numerical techniques in 
 control and other applications,) in general these methods lend to 
 many algorithmic advantages, and the numerical stability which 
 results is important in control problems which can be ill- 
 conditioned. 
 
 The next section discusses the updating of solutions of 
 linear equations when only a few elements of several columns or 
 rows are changed. We will denote this method by "solution per- 
 turbation" since the method depends on solutions to a nominal 
 system of equations. This is a bit of a misnomer, since pertur- 
 
bation implies small changes. However here the changes in ele- 
 ments of the matrix need not be small in any way. In many cases 
 this method has a tremendous computational advantage over factor- 
 ization modification but it does have its drawbacks. This method 
 requires that columns of the nominal inverse, (that is of 
 (I - A) before any elements have changed,) to be computed. In 
 the case of modifying an entire column solution perturbation 
 would require that the entire nominal inverse be computed. Since 
 the computation of an inverse requires approximately three times 
 the computation required to factorize a matrix, this method is 
 not appropriate if many elements in a column are to be modified. 
 Also the method has a disadvantage in that several successive up- 
 dates of the solution due to new perturbations in the matrix ele- 
 ments are difficult to perform since there is no efficient way to 
 compute the updated inverses. 
 
 In Chapter 3 a nonlinear 10 model is proposed which can be 
 solved quite efficiently with the eclectic use the algorithms 
 described in the preceding chapters. This approach largely elim- 
 inates the linearity contraint of the standard 10 model while the 
 extra cost of solving the nonlinear model is quite small. Using 
 solution perturbation the nonlinear equation that must be solved 
 iteratively is reduced to smaller order. Once the residuals are 
 sufficiently small factorization modification will be utilized to 
 find the complete solution. These algorithms have been coded on 
 a minicomputer system, and the computation of solutions of 
 350 -order 10 equations is quite feasible on such a computer 
 system when these algorithms are used. 
 
CHAPTER 1 
 FACTORIZATION MODIFICATIONS 
 SECTION 1.1 
 
 MATRIX DECOMPOSITION BY ELEMENTARY TRANSFORMATIONS, AND 
 SIMPLIFICATIONS THEREOF ON 10 MATRICES. 
 
 This section describes the LU decomposition algorithm. 
 The presentation will also include a discussion of properties of 
 10 matrices that will simplify the algorithms to be described in 
 the next section. 
 
 As with all decomposition algorithms, in the decomposition 
 of the matrix B, (in our case B= I- A, where I is the identity 
 matrix, and A is the technical coefficient matrix defined in the 
 introductory section, thus b^. < and b ii = 1 - a ii > 0,) a 
 transformation is determined such that 
 
 QB= U (1.1.1) 
 
 where U is a upper triangular matrix. To solve the system of 
 equations 
 
 Bx = y (1.1.2) 
 
 we premulitply both sides of (1.1.2) by Q, so that 
 
 QBx = Ux = Qy (1.1.3) 
 
 and the equation Ux = Qy can be solved by back substitution. In 
 the pure form of LU decomposition Q is the composition of only 
 elementary transformations. An elementary transformation has a 
 
8 
 
 simple matrix representation. Its diagonal elements are all one 
 and it has only one nonzero off diagonal element which we shall 
 call the multiplier element. If M i - is an elementary transforma- 
 tion then it looks like: 
 
 m 
 
 ID 
 
 (1.1-4) 
 
 To simplify the notation, the subscripts below a symbol 
 representing an elementary transformation will denote the posi- 
 tion of the multiplier element in the elementary transformation 
 matrix. Thus m. . is in the position (i,j) of the matrix. It is 
 easily verified the the inverse of an elementary transformation 
 is merely the elementary transformation with the nonzero off di- 
 agonal element of opposite sign. The direct multiplication of 
 elementary matrices, 
 
 N. = to .m_ t _: M^, -iM-,, • 
 
 (1.1.5) 
 
looks like: 
 
 j + 2,j 
 
 m 
 
 no 
 
 (1.1.6) 
 
 The LU decomposition algorithm determines the transforma- 
 
 tion 
 
 Q = N n-l N n-2 
 
 N 2 N 1 
 
 (1.1.7) 
 
 If for each j we define 
 
 N j (E j )= N j (N j . 1 N j _ 2 ...N 1 B) , 
 
 (1.1.8) 
 
 then the algorithm determines N. such that it zeroes the elements 
 of the j column of B. below the diagonal element. Thus 
 
 m 
 
 13 
 
 13 
 
 33 
 
 for j + Ki<n 
 
 (1.1.9) 
 
10 
 
 where the n^j's are elements of (1.1.6), B^ = (E^), and d".^ is 
 called a pivot element. 
 
 A reader familiar with decomposition by elementary 
 
 4 
 transformations probably notes that row permutation has been om- 
 itted. The reason for this will be clear shortly. 
 
 Proposition ; If all the column sums of A are less than unity then 
 the pivot elements are nonzero. 
 
 Proof : Clearly if the column sums of A are less than unity then 
 the same must hold for the i leading principal minor of A. 
 Thus by Bellman [1, pp.296, theorem 6] the principal minors of 
 (I-A) are nonsingular. Stewart[10, pp.120, theorem 2.5] has 
 shown that if the principal minors of a matrix are nonzero then 
 if the matrix is decomposed by elementary transformations without 
 pivoting, the pivot elements will be nonzero. 
 
 There is a fundamental result in Input Output analysis 
 which describes the necessary and sufficient condition on A so 
 that the 10 equation admits only positive solutions vectors for 
 positive final demand vectors. This is called the Simon-Hawkins 
 condition [9] . The condition is that all the principal minors of 
 (I-A) must be positive. Let us examine the computation of the 
 solution of 10 equation by factorization so that we can see when 
 
 4 
 
 See Forsythe and Moler[16], or Isaacson and 
 
 Keller [15] on introductory material on decomposition 
 
 methods. 
 
11 
 
 it is impossible for a solution to be a nonpositive vector when 
 the final demand vector is positive. After applying the 
 transformation Q to the final demand vector, since all the multi- 
 plier elements of Q are nonnegative, each element of Qy must be 
 greater than or equal to the corresponding element of y. In the 
 solution of Ux=Qy, since the off diagonal elements are all nonpo- 
 sitive, (thus in the back substitution all sums are on numbers of 
 nonnegative sign,) the only way that a negative element can occur 
 in the solution is if a diagonal element of U is negative. Thus 
 if all the diagonal elements of U are positive, (which implies 
 that all of the multiplier elements of the elementary transforma- 
 tions must be positive,) then it is impossible for the factorized 
 equation to admit a nonpositive solution vector for a positive 
 final demand vector. This is exactly how the Simon-Hawkins con- 
 dition is checked: the determinant of the k principal minor of 
 (I - A) is merely the product of the k topmost elements of the 
 diagonal of U since det|Q|=l. 
 
 In the definition of the transformation Q one order of ap- 
 plying the elementary transformations was given. Since elementa- 
 ry transformations are nearly identity matrices one would expect 
 that these transformations will commute in certain cases, (clear- 
 ly any permutation of the elementary transformations in N. will 
 yield the same transformation.) The elementary transformations 
 can be commuted just as long as no element becomes nonzero that 
 was intentionally zeroed by the application of a previous elemen- 
 tary transformation. Thus: 
 
12 
 
 Q * * 1 nn-l--- M n2"- M 32 M nl--- M 21 
 
 (1.1.10) 
 
 = M nn-1"- M 43^42 K 41 M 32 M 31 M 21 
 
 (1.1.11) 
 
 If we apply Q by (1.1.10) then after the first n-1 elementary 
 transformations are applied, the partially decomposed matrix has 
 the structure: 
 
 X X X X X X X 
 
 X X X X X X 
 
 X X X X X X 
 
 
 
 
 
 u 
 
 (1.1.12) 
 
 After the application of the next n-2 elementary transformations 
 the matrix has the structure: 
 
 
 
 
 
 XXX 
 
 X X X X X X 
 
 X X X X X 
 
 X X X X X 
 
 X X X X X 
 
 X X X X X 
 
 (1.1.13) 
 
 If we apply Q by (1.1.11) then after the application of the first 
 
13 
 
 elementary transformation the partially decomposed matrix has the 
 structure : 
 
 x x x x x x x 
 
 X X X X X X 
 
 (1.1.14) 
 
 Then after the application of the next two elementary transforma- 
 tions the matrix has the structure: 
 
 X X X X X X X 
 X X X X X X 
 X X X X X 
 
 (1.1.15) 
 
 Another ordering which will reduce page faults in virtual 
 computers or 10 requests in successively reading in and writing 
 out parts of the matrix on machines with little main memory in- 
 volves storing the matrix in blocks of rows. The topmost par- 
 tially decomposed block is completely decomposed. Then using this 
 completely decomposed block the remaining blocks are partially 
 
14 
 
 decomposed using this block, that is, all the transformations 
 that need this completely decomposed block will be applied to the 
 lower blocks. Thus this block will no longer be needed. This 
 will reduce page faults or 10 requests by a factor equal to the 
 number of rows that are stored in each block. All these methods 
 are equivalent not only mathematically but numerically, that is, 
 if one is careful about the ordering of the computations, all of 
 the methods will achieve identical matrices in terms of the bit 
 patterns. 
 
 The name of the method, LU decomposition, refers to the 
 definition 
 
 LU = B (1.1.16) 
 
 where L is lower triangular. By inspection 
 
 L = Q" 1 . (1.1.17) 
 
 This is easily shown by the direct multiplication of the inverses 
 
 of the elementary transformations. 
 
 3 
 
 The algorithm performs approximately n /3 multiplications 
 
 and requires only the original matrix as work space. The algo- 
 rithm is the fastest and most compact of all decomposition algo- 
 rithms, Moler [ 11] . 
 
15 
 
 SECTION 1.2 
 
 MULTIPLE ROW AND COLUMN MODIFICATIONS 
 
 In this section an algorithm will pp discussed which will 
 allow the modification of the factorization of a matrix when cer- 
 tain preselected rows and columns of the matrix are changed. A 
 significant advantage of this algorithm is that each modification 
 of the factorization adds no complexity to solving the new system 
 of equations, since the modified factorization is solved in the 
 same manner as the original system. Also the modified factoriza- 
 tion is identical to the factorization resulting from completely 
 factorizing the modified matrix. As we shall see the computation 
 of the modified factorization and the computation of solutions to 
 the modified system of equations can be performed simultaneously. 
 This is useful when the computations are performed on a minicom- 
 puter. Experience has shown that the cost of reading in the ma- 
 trices is the most significant cost in the solution of a large 
 factorized system of equations. Thus streamlining the algorithms 
 with respect to the number of times that the matrices must be 
 read in from secondary storage is worthwhile. In Aoki[4] and 
 Householder [5] , LU decomposition is introduced in a manner that 
 would lend itself to factorization modification methods, though 
 there was no intention to discuss the subject. 
 
 The algorithm for factorization modification greatly sim- 
 plifies if the rows to be modified are at the bottom of the ma- 
 trix, and the columns are on the right hand border. Thus before 
 performing the factorization we exchange the rows to be modified 
 with the bottom rows of the matrix, and exchange the columns to 
 
16 
 
 be modified with the columns on the right hand border. Let P, 
 P~ «P , be the transformation that exchanges the rows and 
 columns. Clearly the matrix representation of P has n^-n zero 
 elements, the other n elements being unity. Rows of the 
 transformation corresponding to coordinates that are not permuted 
 have 1 on the diagonal. If coordinate i is to be exchanged with 
 the j coordinate, then the elements in positions (i,j) and 
 (j,i) would be one. Therefore for all i,j p. .«p.. or P*P t . Also 
 it is easily verified that PP=I. 
 
 For example let there be two industries for which we 
 desire column or row modifications, their positions being say 1 
 and 3. We wish to construct a transformation with the properties 
 above which permutes these industries so that their rows are at 
 the bottom and the columns are on the right hand border. P would 
 then look like: 
 
 
 
 
 
 
 
 
 
 1 
 
 
 
 
 
 1 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 1 
 
 
 
 
 
 
 
 1 
 
 
 
 
 
 1 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 1 
 
 
 
 
 
 
 
 (1.2.1) 
 
 The above example verifies that P is syntpetric. The permuted 
 form of B will be denoted by 
 
 B ■ P BP - PBP . 
 
 (1.2.2) 
 
17 
 
 Let Q be the composition of M.,'s for l<j<n-k, j+l<i<n. 
 Q will partially decompose the matrix B. The elementary 
 transformations are to applied as not to destroy zeros previously 
 introduced. The partial decomposition will be denoted by 
 
 Q p B * . 
 
 (1.2.3) 
 
 In our 6x6 example above, pre- and post-multiplying the original 
 matrix A by P will permute the rows and columns 1 and 3 to the 
 bottom and right border of the matrix. Note that P(I - A) P ■ 
 (I - PAP) . After Q is applied to our permuted (I-A) matrix, the 
 resultant matrix has the structure: 
 
 
 
 
 
 X 
 
 X 
 
 X 
 
 X 
 
 
 
 
 
 
 
 X 
 
 X 
 
 X 
 
 
 
 
 
 
 
 
 
 X 
 
 X 
 
 
 
 
 
 
 
 
 
 X 
 
 X 
 
 (1.2.4) 
 
 At this point any of the last k columns of B can be modi- 
 fied. Let us assume that we want to replace a column of B which 
 corresponds to one of the last k columns of B by a column vector 
 g. is updated merely by replacing the corresponding column of 
 u* by QpPg. Changing more such columns requires the identical pro- 
 cedure for each column. The transformation of g and the right 
 
18 
 
 hand side vector by Q p can be performed simultaneously. This 
 feature in very useful when the computation is performed on a 
 minicomputer system. 
 
 Also any of the last k rows of B can be modified. If a 
 row of B which corresponds to one of the last k rows of B is to 
 be replaced by the row vector h , then the corresponding elemen- 
 tary transformations in Q which zeroed out the first n-k ele- 
 ments of corresponding row of are recomputed. Suppose that in 
 our 6x6 example we desire to replace the first row of B by h . 
 This would correspond to replacing the 5 row of B by (Ph) . 
 Thus we recompute the elementary transformations M_, through Mc d . 
 Before these elementary transformations are applied, the modified 
 has the structure: 
 
 r 
 
 X 
 
 X 
 
 X 
 
 X 
 
 X 
 
 X 
 
 
 
 X 
 
 X 
 
 X 
 
 X 
 
 X 
 
 
 
 
 
 X 
 
 X 
 
 X 
 
 X 
 
 
 
 
 
 
 
 X 
 
 X 
 
 X 
 
 X X X X X X 
 
 X X 
 
 (1.2.5) 
 
 Note that this procedure does not really violate the rule that no 
 element that has been intentionally zeroed be set nonzero by a 
 elementary transformation applied later. Since we are recomputing 
 all the elementary transformations that zeroed out the elements 
 in the row that is to be modified, it is as though these elemen- 
 tary transformations were never applied before the modification 
 
19 
 
 was made. 
 
 For any number ot row and/or column modifications, the 
 modified factorization will be identical to the factorized matrix 
 which would result from partially decomposing B, where B is the 
 matrix with the row and/or column modifications. In fact when 
 the original rows and columns are replaced in the factorization, 
 the new file is exactly bit comparable to the original. To 
 achieve this the reader is warned that care must be taken in the 
 ordering of operations in the original decomposition and the rows 
 and columns modification algorithms. 
 
 To solve the original or a modified system of equations, 
 the rest of the elementary transformations are applied that will 
 reduce to a triangular matrix. This second composition of ele- 
 mentary transformations, (the M^s' where n-k+l<i<n-l and 
 i+l<j<n,) are applied in an order that will not make any element 
 nonzero which was intentionally set to zero by a previously ap- 
 plied elementary transformation. Denoting this transformation by 
 Q and letting TJ be the completely decomposed matrix, the system 
 is solved as follows. 
 
 Bx = y (1.2.6) 
 
 PBPPx = BPx ■ y (1.2.7) 
 
 Q c Q p BPx = Q c UPy = UPx = O c Q p Px (1.2.8) 
 
 Px = «- 1 Q c Q pPy (1.2.9) 
 
20 
 
 x « PU" 1 Q c O p Py (1.2.10) 
 
 Here we have used the fact that PP = I and that P fc * P. 
 The transformation U" 1 is performed by backsubstitution on U. 
 Note that if row and column modifications as above have been per- 
 formed, the algorithm to solve the system does not change. The 
 relevant multiplier elements in the elementary transformations of 
 Q p are changed only if row modifications were performed, while 
 for column modifications only the relevant columns of have dif- 
 ferent elements. If the decomposition is terminated so that the 
 
 lower right hand kxk submatrix remains unfactorized, the modifi- 
 
 1 2 
 
 cation of i rows and j columns requires no more than ^(i+j)n 
 
 multiplications, while the solution of the new system requires 
 
 2 11 
 n + 4k . If the factorization of the kxk submatrix is stored, 
 
 the subsequent solutions require n multiplications. 
 
 To give the reader a sense on how large k can be such that 
 
 the completion of the factorization for each new solution becomes 
 
 comparatively expensive let us set the number of multiplications 
 
 required to solve the system equal to the number required to com- 
 
 3 
 plete the factorization, that is n = ik or k=(3n) For n=100, 
 
 k=45, and for n=400, k=113, and for n=1000, k=208. 
 
21 
 
 CHAPTER 2 
 A SENSITIVITY TRANSFORMATION FOR STUDYING THE PERTURBATIONS OF 
 SOLUTIONS DUE TO THE PERTURBATION OF THE SYSTEM 
 SECTION 2.1 
 
 PERTURBATION OF A SINGLE ELEMENT 
 
 The first example to be investigated is the change of the 
 solution due to the change of one element in the system of equa- 
 tions. This result was first deduced by Shermam and Morrison[8]. 
 Though the result here is the same, the procedure and representa- 
 tion are different. At the end of this section we will point out 
 how this is so. 
 
 Let B be a nxn matrix that is identical to the matrix B 
 except for the element (i,j) which is represented as: 
 
 B ij ■ b ij + 6b ij • (2.1.1) 
 
 Let B be decomposed by a decomposition technique. Then the n 
 columns of the inverse of B can be found by solving the decom- 
 posed system for the appropriate unit vector as a right hand 
 side. Suppose that x is the solution of 
 
 Bx = y (2.1.2) 
 
 and x is the solution of 
 
 Bx = (B+6B)x = y . (2.1.3) 
 
22 
 
 Premultiplying the above equation by B we have 
 
 (I-fB" 1 6B)x « B -1 y * x 
 
 (2.1.4) 
 
 where 6b is a square matrix of order n in which only one element 
 is nonzero, that is, 6 b i-i» Written out explicitly (I + B -1 6B) , 
 
 is of the form: 
 
 B ii* b ij 
 
 B j-i.i 6b ij 
 
 l4 »ji*ij 
 
 b j*i.i 6b ij 
 
 B ni 6b ij 
 
 (2.1.5) 
 
23 
 
 where (B. ,) « B . It follows directly that 
 
 kl 
 
 x. * 13 , (2.1.6) 
 
 1 1+ 5. .6b. 
 
 ji 13 
 
 and for Mj 
 
 x. » x k - b ki 6b, -x, » x k - — *i — iiJ- . (2.1.7) 
 K K K1 1D 3 K 1 ♦ Bj^bij 
 
 This result can also be found by representing the per- 
 turbed system by 
 
 (B ♦ 6b)x * (I ♦ 6BB -1 )Bx = y . (2.1.8) 
 
 If we denote Bx*y as the perturbed system, then solving 
 
 (I + 6BB~ 1 )z=y (2.1.9) 
 
 and then 
 
 Bx = z (2.1.10) 
 
 yields the solution to the perturbed system. The transforming 
 matrix (I + 6BB~ ' is of the form: 
 
24 
 
 Csl 
 
 to 
 
 JO 
 
 \o 
 
 
 w 
 
 
 la 
 
 a 
 
 in 
 
 
25 
 -1 
 
 Solving (I ♦ 6BB~ x )z « y, for Mi 
 
 z k y k , (2.1.12) 
 
 while for the i fc element, 
 
 + (1 ♦ *> lj Bji)i 1 + 6b i;j B ji + 1 y i+1 + ... 
 
 ... + <4b i;( B jn y n , (2.1.13) 
 
 or 
 
 y i " 6b ij B jiyi + 6b ij B j 2 y2 + ••• + 6b ij B ji-iyi-i + 
 
 + (1 + ft^BJDii + »b i;j B ji+1 y i+1 + ... 
 
 ... + «b i3 B jn y n . (2.1.14) 
 
 Since 
 
 n 
 
 * 1 * , ij b jlt»* = 6b ij x j ' (2 ' la5) 
 
 (2.1.14) reduces to 
 
26 
 
 y i " 6b ij\i * (1 * 6b ij b ji )z i " 6b ij 5 jiyi • (2.1.16) 
 
 or 
 
 (1 + 6b ii B ii )y i - ^ ii x i 
 z. U_2i — i iJLJ. (2.1.17) 
 
 1 ♦ 6b ijBji 
 
 6b ii x i 
 
 ■ y, - iJ-2 . (2.1.18) 
 
 1 + 6b. .5. 
 
 Therefore z can be written as 
 
 6b ii x i 
 
 z » y - ( i-i-J )e. (2.1.19) 
 
 1 ♦ 6b ij 5 ji 
 
 where e. is the unit vector of the i component. Since 
 x * B~ 1 z, 
 
 * - x - U-l B 1 e i . (2.1.20) 
 
 1 + 6b. .5 
 
 Explicitly for Mj 
 
 th.46b4.1X 
 
 x k * x. - — !Si — LL_1_ (2.1.21) 
 
 1 ♦ 6 bij 5 D] 
 
and for the j tn element, 
 
 27 
 
 X j * 
 
 b ii 6b Ji X 3 
 1 + 6b ij 5 j . 
 
 1 + 6b i -Bji 
 
 (2.1.22) 
 
 which corresponds to (2.1.6) and (2.1.7). 
 
 The next two sections will consider multiple row and 
 column perturbations. We will see that the analysis of such per- 
 turbations will be quite straight forward. For column perturba- 
 tions the perturbed system of equations will be represented as in 
 (2.1.4), while for row peturbations the representation (2.1.8) is 
 appropriate. In Sherman and Morrison [8] it is only shown that 
 (2.1.6), and (2.1.7) are correct, though they do not show how 
 they arrive at the results. The purpose of this presentation is 
 to demonstrate a procedure which can be used to derive the solu- 
 tion for an arbitrary perturbation, rather than propose a solu- 
 tion and demonstrate its validity. 
 
28 
 
 SECTION 2.2 
 
 COLUMN PERTURBATIONS 
 
 In this section the method just described is extended to 
 columnwise perturbations. The first case is for perturbations in 
 a single column with the scaling of the column, while the second 
 is the extension to two columns. 
 
 Without loss of generality let us suppose that the first k 
 elements of the j column are to be perturbed and that a con- 
 stant multiple of the column is to be added to the entire column, 
 that is, 6b is an nxn matrix which has all zeroes for its ele- 
 ments except for the j column, where the j column is: 
 
 J6b lj +cb ij 6b 2j +cb 2j ... *b kj +cb kj cb k+lj ... cb^l*. 
 
 (2.2.1) 
 
 Again the matrix (I + B~ 1 6B) has the form of (2.1.5) except the 
 
 j column is replaced by: 
 
 where 
 
 J**! * 2 '•• *j-l < 1+ *j +c > *j+l ## ' ^nT (2.2.2) 
 
 '■ = piiV^pj ■ < 2 - 2 - 3 > 
 
 and (B..)* B~ . The solution is found oy inspection as in 
 
29 
 Chapter 2.1: 
 
 and for i;*j # 
 
 *i " x i " *i k j (2.2.5) 
 
 where x is the solution of (B + 6B)x=y, and x is the solution of 
 the unperturbed system Bx*y. 
 
 Now suppose that two columns are to be perturbed , say 
 columns j. and j 2 . For the simplicity of demonstration let us 
 suppose that only the first k 1 and k 2 elements of each respective 
 columns are perturbed. So 6b has only 2 nonzero columns, and the 
 j 1 tn is of the form 
 
 |6b 1H 6b 9 _i ... 6b k • ... | fc (2.2.6) 
 I 1J 1 ^1 K lJl I 
 
 and the j 2 column is: 
 
 K ! t 
 
 |6b,. 6b 9 . ... 6b k . ... r . (2.2.7) 
 
 Now two columns of (I + B~ 1 6B) are not unit vectors, namely 
 columns j. and j 2 . The ji tn column will be denoted by 
 
!*1 *2 
 
 J1-1 
 
 D l 
 
 j,-n # " 
 
 *n| 
 
 30 
 (2.2.8) 
 
 where 
 
 *i ■ p l^9 6b Ph ' 
 
 (2.2.9) 
 
 th 
 
 and the j 2 column as, 
 
 \#1 2 ••* ^ 
 
 j r l (1+ V 'j^l — nj 
 
 (2.2.10) 
 
 where 
 
 0. * 2 5. 6b 
 1 p-1 1P P3 2 
 
 (2.2.11) 
 
 Again by inspection the j* and J 2 th elements of x are 
 
 found by solving 
 
 *| 1+0-i 
 
 (2.2.12) 
 
 Then for Mj 1# j 2 : 
 
31 
 
 *k * x k " V^ " *k 8 j 2 • (2.2.13) 
 
 If all the elements of a column are to be perturbed then 
 it is necessary to have the entire inverse computed. In such a 
 case it would be wiser to use the factorization modification al- 
 gorithm. In many cases thouqh, especially in Input Output 
 analysis, only a few elements are modified, and the rest are just 
 scaled. In such cases the above method becomes very appealing 
 since the simpler closed form equations provide greater insight 
 into the system. 
 
32 
 
 SECTION 2.3 
 
 JKOW PERTURBATIONS 
 
 This section nearly duplicates the presentation of the 
 previous section except that the perturbation vectors are rows 
 instead of columns. 
 
 Without loss of qenerality let us suppose that the first k 
 elements of the j row are to be perturbed and that a constant 
 multiple of the row is to be added to the entire row. Then 6b is 
 an nxn matrix which has all zeroes for its elements except for 
 the j row. The j row is: 
 
 6b jl* cb ji 6b j2 +cb j2 ••' 6b jk +cb jk cb jk+l ••• cb jn 
 
 (2.3.1) 
 
 i-l. , 
 
 Again the matrix (I ♦ 6BB ) is of the same form as (2.1.11) of 
 Chapter 2.1 , but the j row is replaced by: 
 
 where 
 
 I*! * 2 •"• *j-l (1+ *j +c) *j+l" - *ni ' (2.3.2) 
 
 '" = pli 6b 3PV ' (2 - 3 ' 3) 
 
 and (15^)* B" 1 . By inspection the solution of ( I +6BB~ 1 )z= y 
 is, for ij*j: 
 
33 
 z i ■ Yi (2.3.4) 
 
 while 
 
 ^V. + <* + *) + c » z j ■ Vj < 2 -3-5) 
 
 ^.-•A + ^Vj - ry, + (1 + ^ + c)z j (2.3.6) 
 
 pi^jp"? " Vj + (1 + »*j + c,2 j ,2 - 3 - 7) 
 
 and therefore 
 
 ( J + >V y J • D ?/ b 3P X P 
 
 z j " (1 * TC H) ,2 ' 3 - 8) 
 
 where x is the solution of the unperturbed system Bx=y. If we 
 denote z . - y^ as t then 
 
 z = y + te. (2.3.9) 
 
 where e. is the unit vector of the j component. Denoting the 
 solution of the perturbed system as x , as in Chapter 2.1, 
 
34 
 x « B" 1 z * B -1 (y ♦ te.) * x + tB -1 e. . (2.3.10) 
 
 Now suppose that two rows are to be perturbed, say rows j, 
 and j 2 . Again purely for the simplicity of demonstration suppose 
 that only the first k. and k 2 elements ot each respective rows 
 are perturbed. So 6b has only 2 nonzero rows, the ji of which 
 is of the form: 
 
 |6b. , 6b. ... 6b. . ... (2.3.11) 
 
 I *l l 3 1 Z D 1 K 1 I 
 
 and the j 2 row is: 
 
 |6b. , 6b, ... 6b, „ ... (2.3.12) 
 
 3 2 l J2 Z 3 2 K 2 I 
 
 Now two rows of ( I + 6BB ' are not unit vectors, namely rows j, 
 and j 2# The ji row will be denoted by 
 
 I I 
 
 1 ^2 
 
 r l (1+ >V ^j 1+ l-'- *n 
 
 (2.3.13) 
 
 where 
 
 k i 
 
 *•» = pii' b i ,p B pi ' ,2 - 3 - 14) 
 
th 
 
 and the j 2 row as » 
 
 where 
 
 Sol' 
 
 ,-1 
 
 ving ( I +6bb ) z ■ y, for i^j , j 
 
 while for i= Ji»Jo : 
 
 where 
 
 and 
 
 35 
 
 *1 *2 •'• ^j r l Cl^j ) jl+1 
 
 • • • 
 
 n 
 
 (2.3.15) 
 
 0, - 2 6b. T5 . . 
 1 p=l D 2 p pi 
 
 (2.3.16) 
 
 z . ■ y . 
 3 Y 3 
 
 (2.3.17) 
 
 I 14*, >*, I |z. 
 
 I 
 
 2 II J l| 
 
 0^ 1+0^ !z 
 
 D 2 M 3? 
 I I 2 
 
 (2.3.18) 
 
 y+ - (l ♦ A* )y, + r*, y-i - 2 6b. x 
 
 3i '3i 3o'3 
 
 1 J l J 2 J 2 p=l J l 
 
 JiP P 
 
 (2.3.19) 
 
36 
 
 3 2 
 
 Ya - (1 + *a )y i + 0t Vt - 2 <*b. x D . (2.3.20) 
 J 2 J 2 D 2 D l 3 1 p=l 3 2 p p 
 
 Redefining t to be z. - y. and defininq u to be z^ - y+ , then 
 
 3 1 3 1 D 2 D 2 
 
 the perturbed solution is: 
 
 x = B" 1 z = B" 1 y + tB" 1 e. + uB* 1 e. . (2.3.21) 
 
 D l J 2 
 
 If an entire row is perturbed then only its associated 
 column of the inverse is needed in the computation of the per- 
 turbed solution vector. The computation requires 2n+l multipli- 
 cations, (column perturbation would require about n .) 
 
 The above presentation lends particular insight into the 
 perturbation of row elements in linear equations. When several 
 rows are perturbed then the perturbation of the solution is a 
 linear combination of the associated columns of the nominal in- 
 verse. Therefore we can make a rough quess on whether solutions 
 are sensitive to perturbations of several rows of the 10 matrix 
 by inspecting the associated columns of the n^ninal inverse. 
 
37 
 
 CHAPTER 3 
 A NONLINEAR INPUT OUTPUT MODEL 
 
 As we have seen both factorization modification and solu- 
 tion perturbation have comparative advantages in the solution of 
 modified linear equations. Factorization modification has the ad- 
 vantage of simplicity of data storage, that is, when a row or 
 column is changed, then only that row or column of the modified 
 factorization is changed. Also the repetitive modification of 
 
 the linear equations is very stable numerically since the refac- 
 
 5 
 torized matrix is machine identical to the factorized matrix 
 
 when the modified system of linear equations is completely decom- 
 posed from scratch. The disadvantage of factorization modifica- 
 
 2 2 
 
 tion is that n multiplications must be performed, and n matrix 
 
 elements must be read in from secondary storage. In some cases 
 solution perturbation requires much less computation than factor- 
 ization modification to achieve the same result. The disadvan- 
 tage of solution perturbation is that if many elements in dif- 
 ferent rows are to be changed, as in the case of modifying all 
 the elements in a single column, then much or all of the nominal 
 inverse must be computed. Also repetitively computing new solu- 
 tions due to changes in the elements of the matrix is difficult 
 by this method. In this section we will utilize the comparative 
 advantages of both methods in solving a nonlinear Input Output 
 model . 
 
 5 
 
 "Machine identical" means that the bit patterns of 
 
 the two matrices are identical. 
 
38 
 
 There has been very little work presented in the litera- 
 ture about nonlinear 10 models, most of it only in the way of ex- 
 istence proofs, Sandberq[6) and Lahiri[?J, and V3ry little empir- 
 ical application. This may be due to the difficulty of obtaining 
 data on the nonlinear ity of the transactions, for just obtaining 
 the transactions matrix data is a difficult task in itself. But 
 the nonlinear ity of the transaction elements may not be due so 
 much to the nonlinearity of the input requirements per unit out- 
 put for each firm. Rather the nonlinearity may be due to the 
 average input requirements per unit output, that is the a^'s, 
 shift from the production techniques of firms that cannot in- 
 crease output towards the production techniques of those firms 
 that can. A good example of this type of nonlinearity is the 
 transactions to the oil industry if oil demand would change. If 
 the demand for oil by consumers and industry would drop drasti- 
 cally then the transactions for drilling equipment and real 
 estate would drop disporportionally to the drop in total oil out- 
 put since oil wells would last for a longer period of time. 
 Given more time there would be fewer chances that wells that are 
 drilled end up to be dry, etc. It is likely that at lower oil 
 demand more oil drilling would be internally financed, and thus 
 the amount of interest charges would be proportionally less. 
 
 The transactions matrix is AX where X is a 
 diagonalization of the total demand vector. The 
 transactions matrix is the interindustry data that 
 governments actually collect, the A matrix is then derived 
 from the transactions matrix after the total demand vector 
 is computed by summing the rows of the transaction matrix 
 and adding the final demand vector to the resultant vector. 
 
39 
 
 Without scarcity driving up prices firms are less willing to take 
 chances, they may curtail research and development activity. 
 Thus the composition of value added changes. As we shall see in 
 the example below, if the demand for oil would rise the a. . # s for 
 the oil industry would shift toward the production techniques of 
 the exploratory firms. New production techniques such as extrac- 
 tion of oil from shale would become more profitable, and average 
 production shifts towards this technique. By classification of 
 firms into two groups, those that can and cannot increase output, 
 we have an elementary procedure by which we can represent the 
 nonlinearity of transactions. 
 
 Let j be the index corresponding to the oil industry. 
 Suppose that oil total output is represented by 
 
 X j * Xn j + X °j ' (4.1.1) 
 
 where x . , x n • , x°^ are total demand for oil, total demand for oil 
 supplied by new oil, and total demand for oil supplied by old oil 
 respectively. We assume that only the producers of new oil can 
 expand the output of oil to a level greater than x .. The tran- 
 saction from the i industry to the oil industry is 
 
 t tj = a ijXj = a° ij x° j + a n ij(Xj - x^) . (4.1.2) 
 
 Therefore 
 
40 
 
 / o .n , v o 
 t ii < a ii ii ,x ii n 
 13* r x^ x j r ij H.i.JJ 
 
 for x.^x ^ . (4.1.3) exemplifies the shift in the technical coef- 
 ficients mentioned in the preceding discussion. 
 
 Though oversimplified, (but the linear 10 model is even 
 more simplified,) this example does illustrate that empirical 
 work in this area may be more practical than was previously 
 thought. It can also open up the possibility of treating invest- 
 ment as a direct input to production instead of treating it as 
 final demand. Investment cannot be realistically treated as a 
 linear function of total output since if existing capital is used 
 at full capacity the only way output can increase is by invest- 
 ment. Thus investment would make up an considerable portion of 
 costs. Alternatively if capital is not being used at full capa- 
 city then there will be little investment till excess capacity 
 has depreciated away. Observation of plant capacity levels, and 
 the catagor ization of firms in a industry may give us the essen- 
 tial information needed to produce a viable nonlinear 10 model. 
 
 The nonlinear 10 model in the following exposition has an 
 restriction in the form of the nonlinear ities. It is not the most 
 general model that could be presented, though the formulation 
 does encompass the situations that we would most likely encounter 
 in the real world. We assume that the nonlinearities of the 
 tranasctions are only a tunction of the buying industries, that 
 is 
 
41 
 t ij * t iJ (X J ) " (4.1.4) 
 
 This is a crucial assumption in this algorithm for it tremendous- 
 ly reduces the order of the computation. Though this is a res- 
 triction to the analysis, it is difficult to visualize an example 
 where a transaction element is a direct function of total demand 
 other than the buying industry. Even in the literature though 
 this assumption is seldom made or utilized, the examples present- 
 ed usually are of this form. Also we have assumed a functional 
 form for the nonlinear it ies. Though this assumption is not cru- 
 cial to the exposition, and it can easily be relaxed, we shall 
 see that this representation is computationally advantageous. 
 
 A truncated power series will be used as the functional 
 form of the nonlinearities. The use of only three terms is only 
 for the simplicity of demonstration. 
 
 x. x. 
 
 (4.1.5) 
 
 where 
 
 i ■ «ij + Pij + 'ij 
 
 (4.1.6) 
 
 Note that 3i-j =a i-i (x^) . We denote the matrix A" and vectors x and 
 7 as the nominal technical coefficients matrix, »nd nominal total 
 and final demand vectors respectively. If c<. .=1, and Pi-\~Yi-h s ^ 
 
42 
 
 for some a j-i(x.i)» then the function is constant, as in the linear 
 10 model. Varying the values of p i - and Y i - from zero introduces 
 nonlinear ity in the technical coefficient. But as long as con- 
 straint (4.1.6) holds then the nominal output vector x is still 
 the solution of the nonlinear 10 model when the final demand vec- 
 tor is y. 
 
 To introduce the concept of the model, let us assume that 
 only the j industry has nonconstant technical coefficients. 
 These technical coefficients are only a function of the total 
 output of the j industry, that is, x.. Let y be a vector dif- 
 ferent from y and let x 1 be the solution of 
 
 (I - T^x 1 = y . (4.1.7) 
 
 Clearly there is little chance that x 1 is the solution of 
 the nonlinear equation 
 
 (I - A(x))x = y . (4.1.8) 
 
 Let 
 
 AMx) = A(x) - A . (4.1.9) 
 
 Note that AA(x) is nonzero only for the j column. By solution 
 perturbation 
 
 (I - A(x))x = (I - A)(I - (I - AJ^AAfx))* = y (4.1.10) 
 
43 
 
 or 
 
 (I - (I - A) _1 ^(x))x = x 1 (4.1.11) 
 
 where x is the solution of (4.1.7) 
 
 The j th row of 4.1.11 is 
 
 n 
 
 [1 - I 5,-6 a.-fx,)]*, = x* (4.1.12) 
 
 i = l J x a j J J J 
 
 where 
 
 -1 
 
 (5^) = (I - A)" 1 , (4.1.13) 
 
 and 
 
 (^a ij (x j )) =AA(x) . (4.1.14) 
 
 Since the x. is the only unknown in (4.1.12) we need only solve 
 this scalar nonlinear equation for x.. If the nonlinearities are 
 represented as a truncated power series as in (4.1.5) , then 
 (4.1.12) reduces to 
 
X • X • 
 
 [1 - en-*) - 1) - a((-J-) 2 - l)]x = x* 
 x. I 
 
 44 
 
 (4.1.15) 
 
 where 
 
 n 
 
 * s 2 ^iiaiiPii r and (4.1.16) 
 
 e = 2 B.x. y.. . (4.1.17) 
 
 i=l D1 13 J 
 
 Then (4.1.15) can be solved by a numerical method such as 
 Newton's method. 
 
 The advantaqe of the functional representation (4.1.5) is 
 that the effects of the nonlinear ities of all the elements of the 
 j column of A(x) are expressed in cr and 0. If more terms are 
 added to (4.1.5) then only similar terms are then added to 
 (4.1.15). Note that only the j row of the inverse needs to be 
 computed to solve for x . , but once we have solved for x . we now 
 know all the values of the technical coefficients in the j 
 column. Using factorization modification we can compute the en- 
 tire solution without the computing the inverse, as ve would have 
 had to do if we use solution perturbation only. During the fac- 
 torization modification we can simultaneously compute the solu- 
 tion x, thus saving an extra pass through the matrix. ^ 
 
 3 
 
 This is only important of course if the computer 
 
 installation does not have enough main memory to hold the 
 
 entire matrix. 
 
45 
 
 The extension to case of nonlinearities in more than one 
 column of A is quite straightforward. If k columns of A have non- 
 linearities, then the k corresponding rows of (I - "R)~ must be 
 computed and a k order nonlinear matrix equation corresponding 
 to (4.1.15) is solved. 
 
 This algorithm is a tremendous savings over the usual 
 Newton's method solution of this problem, which would be 
 
 x k+1 = x k - (I - ^_T)" 1 [x k - A(x k )x k - y] (4.1.18) 
 
 k t" h 
 
 where x is the approximate solution at the k iteration, and T 
 is a vector valued function whose j element is 
 
 t i « 2 a..(x.)x. . (4.1.19) 
 
 For each iteration (4.1.18) requires approximately 
 
 3 2 
 
 n /3 + 2n + ikn multiplications where n is the order, i+1 is the 
 
 number of terms in the truncated power series, and k is the 
 number of nonlinear columns. The method proposed here requires 
 only about k /3 + (i+l)k multiplications. If n=100, k=5, and 
 i=3, then an iteration of (4.1.18) requires about 353,000 multi- 
 plications, while the method presented here requires about 162 
 multiplications. Thus the proposed method has a computational 
 savings by a factor of 2,180! 
 
 Most importantly this method gives us a very simple means 
 to compute the "marginal inverse", that is compute the incremen- 
 
46 
 
 tal total output due to an increment in final demand. The margi- 
 nal inverse is the inverse of the Jacobian matrix evaluated at a 
 solution. Sandberq [6, Theorem 1] has shown that the inverse of 
 the Jacobian evaluated at a solution will approximate the pertur- 
 bation of solutions around the solution due to perturbations of 
 final demand. If the vector c is the perturbation in final 
 demand, and the vector x^ is the actual perturbation of total 
 demand then 
 
 x p = J" 1 c + A»(c) (4.1.20) 
 
 where 
 
 1 |Aj£j ! ! -> as ||c|| -> . (4.1.21) 
 
 The norm operator is the Euclidean. In our example in which only 
 column j is nonlinear, if (J ik ) = J, then: 
 
 j ik = -a ik , for k*i,j (4.1.22) 
 
 j Rk = 1 - a Rk ^ for Mj (4.1.23) 
 
 and 
 
 Jij = " a ij " ( x j)^5T7 a(x j ) for i ^' (4.1.24) 
 
47 
 
 j jj = * ~ *jj " (X j ) ^T7 a(X j ) ' (4.1.25) 
 
 Since the Jacobian matrix is identical to the 10 matrix 
 except tor the j column we can use factorization modification 
 to replace the j column of the nominal 10 matrix by the the as- 
 sociated column of the Jacobian matrix. In doing so we have ef- 
 fectively computed the factorization of the Jacobian. 
 
 At this point it is worthwhile to digress to a fundamental 
 result of 10 analysis. In the linear 10 model price is computed 
 by the identity 
 
 (I - A fc )p = v or p = A fc p + v . (4.1.26) 
 
 The j row of (4.1.26) reads: the price of good j, that 
 is, p. equals the average unit costs of inputs, that is, (A p) . 
 plus the cost of value added, that is, v.. In this equation pro- 
 fit rates must be assumed, and the costs are average, not margi- 
 ns 
 
 nal. Assuming nonincreasing returns suppose that we are given a 
 final demand vector y and the corresponding total demand vector x 
 which is the solution of the nonlinear 10 model. Then we can 
 find a price vector p such that it is profit maximizing for each 
 industry to produce the total demand vector x when the final 
 
 o 
 
 Concavity of the technical coefficients functions and 
 value added functions would imply nonincreasing returns. 
 
4t 
 
 aemand vector is y. Tne equation wnich computes tnis price level 
 is : 
 
 J P - v 
 
 (4.1.27) 
 
 wnere cue i element ot v" is 
 
 v(x.) 4,x i ) 5 H-v(x 1 ) . 
 
 1 
 
 (4.1.2b) 
 
 we have allowed tor nonlinear ities in value added, and as tor the 
 transaction elements, the nonlinear ities are only a runction ot 
 that particular industry. ol course tnis definition ot value 
 added excludes profits. we can easily verity that the price 
 vector which is the solution ot (4.1.27) is the price vector that 
 makes x the profit maximizing output vector. tor the j 
 industry tne per unit profit is 
 
 n 
 
 2 a. (x • ) - v (x- 
 i=l J J 3 
 
 4.1. 2b) 
 
 'iherelore total profits are: 
 
 n 
 
 " 3 " V P J " i ! 1 a iJ (X 3> - V(X D )] • 
 
 (4.1.3k,) 
 
xaKinq the derivative with respect to x we have 
 
 j 
 
 n 
 
 cJx- 77 j s ^ Pj * .2 1 P i a ij (x j ) + v(x 
 
 J l-l 3 1 
 
 ) J 
 
 (4.1.31) 
 
 or 
 
 'j ' J 1 P i [a ij (x j ) + < x j>c£r*ij<*j>] 
 
 V(Xj) 4 (Xj)^-V(Xj) . 
 
 (4.1.32) 
 
 Notice that (4.1.24) and (4.1.2b) are the general form ot the 
 elements ot the Jacobian matrix J. Also the right hand side of 
 (4.1.32) corresponds to (4.1.2b). by differentiating each of the 
 industry's profits function with respect to that industry's total 
 output we see that (4.1.26) is tne first order condition tor 
 profit maximization. 
 
50 
 
 DISCUSSION AND CONCLUSIONS 
 
 A summary ot the relative advantages and disadvantages of 
 solving modified systems of equations by factorization modifca- 
 tion using elementary transformations and by solution perturba- 
 tion was presented in the beginning ot Chapter 3. The conclusion 
 of the discussion was that neither method had a clearcut advan- 
 tage over the other. When the two methods are used in conjunc- 
 tion with each other, as in the nonlinear 10 model, the resulting 
 algorithm can be very efficient. 
 
 The efficient solution of 10 equations has been largely 
 ignored by economists. This is evident by the fact that up to 
 now only large batch computer systems were used to solve 10 equa- 
 tions. Using the algorithms discussed in this paper a 368 ord- 
 
 Q 
 
 er 10 matrix was factor ized on a minicomputer installation. The 
 solution computation of the solution for an arbitrary final 
 
 demand vector required 8 seconds user time and 12 seconds system 
 
 n ma 10 
 time . 
 
 9 
 
 The computer was a Digital Electronics Corporation 
 
 PDP11/50. The computer was operating un^er the UNIX 
 operating system developed by Bell Laboratories. The 
 algorithms were coded in C, the principle language in which 
 much of the UNIX system was coded. 
 
 10 
 
 The Unix operating system indicates two types of 
 
 program timinqs, one for user time, which is the time 
 
 required to perform the actual computation in the user's 
 
 program area. The other is denoted as system time which is 
 
 the computation required by the operating system to perform 
 
 principly the input output functions. Since a 368 tn order 
 
 matrix requires approximately 1 megabyte of storage, the 
 
 computation of physical block addresses of the data consumed 
 
 the largest portion of the total computation. System time 
 
 would be significantly reduced by performing raw, ie., pure 
 
 direct memory access (DMA) input output. 
 
51 
 
 Double precision arithmetic was used by both the decompo- 
 sition and solution programs. As a test of the numerical stabil- 
 ity of the algorithms on actual data, the Bureau of Economic 
 Analysis (BEA) 368 th order 10 matrix was decomposed, and then a 
 
 system of equations utilizing this matrix was solved. The resi- 
 
 -14 . . 
 
 duals were of the order of 10 when double precision arithmetic 
 
 was used. Thus the use of elementary transformations has not 
 been detrimental to numerical stability. 
 
 As a test of the nonlinear 10 model, nonconstant technical 
 coefficients were introduced into five columns of the 1967 368 
 order BEA 10 matrix. The iterative solution perturbation algo- 
 rithm required 3.3 seconds user and 6.1 seconds system time, 
 while the factorization algorithm which computes the entire solu- 
 tion vector required 31 seconds user and 27 seconds system. 
 Since the solution perturbation algorithm computes the total out- 
 put for the industries which have nonlinear input technical 
 coef f f icients , we can check the computation of the factorization 
 algorithm by comDaring the elements in the solution vector to 
 
 those computed by the solution perturbation algorithm. The 
 
 -14 
 residuals were of the order ot 10 
 
 For a small set of nonlinear industries, the computation 
 involved in the solution of a nonlinear 10 model is about two or 
 three times the computation involved in solving a factorized sys- 
 tem of equations. In terms of computation, there is little that 
 bars the incorporation of the algorithms in empirical 10 
 
52 
 
 research. It may be futile to consider the construction of a em- 
 pirical nonlinear 10 model with all ot the columns being non- 
 linear, at least this is so tor the present. Many times 
 researchers who utilize 10 models tocus most of their attention 
 upon one or two, or a small qroup of industries. The effects of 
 alternative technical coefficients ot a small group of industries 
 are often analyzed. ihe framework presented could easily be 
 molded to such applications. Finally it is suggested that inter- 
 industry data collection should gear some ot its efforts toward 
 the determination ot capacity levels of the individual firms. 
 Doing so will allow the construction of the nonlinear investment 
 functions needed to make the endogenous determination of invest- 
 ment levels in 10 models possible. 
 
53 
 
 REFERENCES 
 
 [I] R. Eellman, I ntr oduct ion to M atrix A nalysis . New York: 
 McGraw-Hill, 1970", 2nd ed. 
 
 [2] A. Sameh and P. Bezdek, "Methods for Increasing the Computa- 
 tional Efficiency of Input-Output and Related Large Scale Matrix 
 Operations," Center for Advanced Computations Document No. 66, 
 Univ. of Illinois, May 1973. 
 
 [3] K. Noh and A. Sameh, "Computational Techniques for Input- 
 Output Econometric Models," Center for Advanced Computations Do- 
 cument No. 134, Univ. of Illinois, Sept 1974. 
 
 [4] M. Aoki, Introduction to Optimization Techniques . New York: 
 Macmillan, 1971. 
 
 [5] A. Householder, The Theory of Matr ices in Numer ical Analysis . 
 New York: Blaisdell, 1964. 
 
 [6] I. Sandberg, "A Nonlinear Multisectored Input-Output Model of 
 a Multisectored Economy," Econometr ica , Vol. 41, No. 6, pp. 
 1167-1182, Nov. 1973. 
 
 [7] S. Lahari, "Input-Output Analysis with Scale-Dependent Coef- 
 ficients," Econometrica, Vol. 44, No. 5, pp 947-961, Sept. 1976. 
 
 [8] J. Sherman and W. Morrison, "Adjustment of an Inverse Matrix 
 Corresponding to a Change in One Element of a Given Matrix," An- 
 nals of Mathematical Statistics, No.l, Vol. 21, pp 124-126, March 
 1950. 
 
 [9] D. Hawkins and H. Simon, "Note: Some Conditions of Ma- 
 croeconomic Stability," Econometrica, Vol. 17, No. 3, pp. 
 245-248, July 1949. 
 
 [10] G. W. Stewart, I ntroduction to Matrix Comput a tion . New 
 York: Academic, 1973. 
 
 [II] C. Moler "Matrix computations with Fortran and Paging," Com- 
 munications of the Association of Computing Machinery, Vol. 15, 
 pp. 268-270, April 1972. 
 
 [12] G. H. Golub, P. E. Gill, w. Murray, and M. A. Saunders, 
 "Methods for Modifying Matrix Factorization," Stanford Univ, Rep. 
 STAN-CS-72-322, 1972. 
 
54 
 
 [13] P. E. Gill and W. Murray, "A Numerically Stable Form of the 
 Simplex Algorithm," Linear Alqebra and Its Applications., Vol. 7, 
 pp. 99-138, 1973. 
 
 [14] G. J. Bierman "Additional Comments on "Multistage Least- 
 Squares Parameter Estimation," IEEE Trans. Automat.. Contr. Vol. 
 AC21, pp. 883-885, Dec. 1976. 
 
 [15] E. Isaacson and H. Keller, Analysis of Numerical Methods . 
 New York: Wiley and Sons, 1966. 
 
 [16] G. Forsythe and C. Moler , Computer Solutions of Linear 
 Algebraic Sys t ems , Englewood Cliffs, N.J.: Prentice Hall, 1967. 
 
 [17] H. Chenary and P. Clark, Inter industry Economics , New York: 
 Wiley, 1959. 
 
UNIVERSITY OF ILLINOIS URB™ 
 
 510.84IL63C C0D1 * 
 
 _3 01l2 0Q72fi4nfiQ