The Method of Lagrange Multipliers

1 The Method of Lagrange Multipliers Phil Lucht Rimrock Digital Technology, Salt Lake City, Utah 84103 last update: June 22, 2015 Maple code is available upon request. Co mments and errata are welcome. The material in this document is copyrighted by the author. The graphics look ratty in Windows Adobe PDF viewers wh en not scaled up, but look just fine in this excellent freeware viewer: http://www.tracker-software.com/pdf-xchange-products-co mparison-chart . The table of contents has live links. Most PDF view ers provide these links as bookmarks on the left. Overview an d Summary ........................................................................................................... .............. 2 1. The R matrix and its Rank................................................................................................... .............. 3 2. The Method of Lagrange Multipliers.......................................................................................... ......6 3. Proof of Theorem 1 .......................................................................................................... ................. 10 (a) Proof for N = 6 and S = 4 .................................................................................................. ............ 10 (b) Proof for general S ≤ N ................................................................................................................. 14 4. The Gradient Interpret ation: Examples 1 and 2............................................................................ 18 4.1 Setup ............................................................................................................................................. 18 4.2 Example 1: Half sphere with one simple constraint..................................................................... 19 4.3 Example 2: Half 4D sphere with two simple constraints............................................................. 24 4.4 The General Case........................................................................................................... ............... 27 5. Example 3: Extremal distan ces between an ellipse and a circle .................................................. 28 6. Example 4: The Boltzmann factor in Statistical Mechanics ......................................................... 35 6.1 Statement of The Boltzmann Extremum Problem ........................................................................ 35 6.2 Solution of The Boltzmann Extremum Problem........................................................................... 37 6.3 More details of the solution ............................................................................................... ........... 40 6.4 A numerical example with plots ............................................................................................. ......43 Appendix A : Matrix rank equals the number of independent columns and rows......................... 49 Appendix B : Determinants.................................................................................................... ............ 53 B.1 Definition of the determinant .............................................................................................. ......... 53 B.2 The permutation group a nd the permutation tensor ε................................................................... 55 B.3 Minors, Cofactors, and the Cofactor Expansions of det(M) ........................................................ 59 B.4 Expressions for the inverse matrix M-1 and Cramer's Rule ......................................................... 61 References.............................................................................................................................................. 64 2 Overview and Summary The Method of Lagrange Multipliers is used to determine extrema of a real function subject to some number of constraint equations. The main purpose of this document is to provide a solid derivation of the method and thus to show why the method works. A sec ondary purpose is to provide some interesting and illustrative examples. Although we discuss the relevant gradient approach in Section 4, it seems more convincing to derive the method using a matrix approach where the candi date solution points of a constrained extremum problem are associated with the points at which a cer tain matrix drops below full rank (an approach used by Buck, see References). Then the gradient stat ement of "why it works" follows at once. Every document on this subject needs a "Theorem 1" and ours is no exception. Fortunately for the reader, there are no other theorems (ignoring the appe ndices), and Theorem 1 is fairly easy to prove. The reader should be at least dimly familiar with the basic facts of linear algebra. Here is a concise summary of the document: Section 1 gives a statement of Theorem 1 in the context of a certain R matrix Section 2 shows how Theorem 1 provides an immediate de rivation of the Method of Lagrange Multipliers. Section 3 then proves Theorem 1. Section 4 gives the gradient interpretation with two visualizable examples. Section 5 presents another example which, though si mple, requires numeric work to solve. Section 6 provides an extended standard example involving statistical mechanics. Appendix A proves for a general matrix that column rank = row rank = rank. Appendix B derives a batch of facts about determinants which are used in Appendix A. References (a very short list) are provided on the last page. The general tone is more that of an engineer or physicist, not that of an abstract mathematician. Maple code is used where it seems useful to implement calculations or display graphs. When an earlier equation is quoted, its equation number is put in italics. 3 1. The R matrix and its Rank We shall operate in N dimensional space EN with coordinates r = (x1, x2 ... xN). We seek an extremum of a real function f( r) subject to S-1 constraints a( r) = 0, b( r) = 0, ... q( r) = 0 where we use q generically to represent the last constraint function. For example, if there are two constraints, then S = 3 and the two constraint equations are a( r) = 0 and b( r) = 0. So here is the problem: find extremum of f( r) subject to constraints a( r)=0, b( r)=0, c(r )=0, .... q(r )=0 . (1.1) 1 2 3 S-1 Construct the following matrix of partial derivatives, where f i(r) ≡ ∂if(r) = ∂ f(r)/∂xi. That is, the subscript on f indicates which of the arguments x 1,x2, x3....xN the partial is with respect to. So here is that matrix, which we shall call R : (1.2) For example, f 2 ≡ ∂f(x1,x2, x3....xN)/∂x2. On the far right we abbreviate each row of functions in bold , indicating a row vector like f . The matrix has S rows and N columns. We shall only be interested in S ≤ N (see comment below). If S = N, the matrix is square, otherwise it has more co lumns than rows, so has a horizontal band shape. The maximum possible rank of the above matrix is S since this is the smaller of the number of rows and columns (see below). Since all the functions are functions of r , we really have R( r). As the point r is varied, the elements of the matrix R( r) generally change. If S = N, the matrix is square and we can talk about the determinant det(R). In this case there is only one SxS "submatrix" and it is the entire matrix R. If S < N, the matrix is wider than it is tall. In this case, we can talk about square submatrices within R which are of shape SxS. For example, for N = 6 and S = 4 we have three constraints and the R matrix is this, (1.3) 4 Here we have outlined two 4x4 submatrices in red and blue. The columns of a submatrix need not be contiguous, so the three green boxes indicate anothe r 4x4 submatrix. In this example, the number of submatrices is given by the binomial factor (6,4) -- pick a committee of 4 columns from 6 candidates. So we have shown only 3 out of (6,4) = 6!/(4!2!) = 15 possible 4x4 submatrices. We can indicate a submatrix using this notation: red = (f 1,f2,f3,f4) blue = (f 2,f3,f4,f5) green = (f 1,f2,f4,f6) (1.4) That is to say, we abbreviate the submatrix by sta ting only the values in the first row of the submatrix. Since our main interest will be whether or not the determinants of these 4x4 submatrices vanish, the ordering of the columns within a submatrix is not important. For the general R matrix shown above with N column s and S rows, there are (N,S) possible submatrices of size SxS. A particular submatrix is then indicated by a sequence of S subscripted f values. Why not S>N? One can regard u = f(x 1,x2, x3....xN) as an N-dimensional surface in EN+1 where the extra dimension is u. For example, for N = 2, u = f(x,y) is a 2-dimensional surface in E3. The constraint equation a(x 1,x2, x3....xN) = 0 is an N-1 dimensional surface in EN. Again for N = 2, a(x,y) = 0 is a 1-dimensional surface (a curve) in E2. If a(x,y) = x2-y2-1, that curve is a circle. A constraint surface can be "extruded" in the u direction of EN to give an N dimensional surface in EN+1. For N=2, if the constraint is a circle , that extruded surface is a cylinder in E3. The intersection of the N-dimensional surface u = f with the N-dimensional extruded surface a = 0 is a surface of dimension N-1. In our N = 2 example, the 2D surface u = f intersects the cylinder in a curve. Each extra constraint lowers the dimensionality of the final intersection surface. If S = 2 (one constraint), the intersection surface was just seen to be of dimension N-1. For S = 3 (two constraints), that surface is down to dimension N-2. The dimension of th e final intersection surface is thus N-S+1. When S = N, the final surface is a 1D curve. If we try to set S > N, the intersection has no degree of freedom at all and one cannot have a continuous extremum situation. An N = 2 example of this paragraph's discu ssion is presented in Section 4 (Example 1). We now quote a few facts from linear algebr a concerning the rank of a matrix: Definitions : A minor of matrix A is the determinant of any square submatrix of A obtained by crossing out some rows and columns The rank r of an m x n matrix A is the dimension of the largest non- vanishing minor within A. [Shilov 1.92 ] Thus, r ≤ min(m,n). (1.5) Definitions : [Shilov 2.21] Consider an m x n matrix A. Let r1, r2, ...rN be an arbitrary selection of N rows (row vectors) from the matrix A. This set of N rows is linearly dependent if one can find a set of constants k i at least one of which is non-zero such that k1 r1 + k2 r2 + .... + k N rN = 0 . ( 1 . 6 ) There could be several such equations, or there could be only one. This generally means that at least one row can be expressed as a linear combination of the other rows. A legal such equation would be r2 = 0 in which case k 2 = 1. One could then say that r2 is a linear combination of the other rows with coefficients all zero. In any event, if one or more rows are all zeros, a set of N rows are linearly dependent. If no such 5 equation (1.6) exists, then the set of N rows is linearly independent . For example, the set of vectors r1 = (1,0) and r 2 = (0,1) is linearly independent. The above discussion also applies to columns. Take "row" → "column" and ri → ci . Facts : The rank r as defined in (1.5) has these prope rties for any m x n matrix A: (1.7) (a) the rank r equals the number of linearly independent columns of A [Shilov 3.12 ] (b) the rank r equals the number of linearly independent rows of A [Shilov 3.13 ] See Appendix A below for a proof of the claims of (1.7). With the above as background, the following theo rem is fundamental to our Lagrange Multiplier presentation: Theorem 1 . If the point r is a solution to the extremum problem (1.1), then rank[R( r)] < S. (1.8) According to the rank definition (1.5), rank(R) < S ⇔ all SxS submatrices must have zero determinant (1.9) If, for example, the red submatrix in Fig 1.3 had a non-zero determinant, then rank(R) = 4 = S. The point of Theorem 1 is that the solution point r of the constrained extremum problem must be a point where the R( r) matrix drops below full rank. We shall prove this theorem in Section 3 by showing that, if r is a solution to the extremum problem, then a ll those SxS submatrices discussed above have zero determinant and thus rank(R) < S. If one is searching for candidate solutions r to the extremum problem, one can restrict one's search to points r for which rank[R( r)] < S. Comment on the R matrix. Imagine a general transformation from x-space x = (x1, x2, ....xN) to u-space with coordinates u = (u1, u2, ... uS). One might write such a transformation as u = F(x) where F:EN→ ES. Although this transformation in general is non-linear, in a tiny neighborhood of point x in x-space (and the corresponding point u in u-space), the transformation will be linear if certain conditions are met. That local linear relation is then described by d u = Rd x where R is a matrix (N co lumns, S rows) whose matrix elements are R ij = ∂ui/∂xj. If we think of f = u 1(x1, x2, ....xN) R 11 = ∂u1/∂x1 = ∂f/∂x1 = f1 R 12 = f2 etc. a = u 2(x1, x2, ....xN) R 21 = ∂u2/∂x1 = ∂a/∂x1 = a1 R 22 = a2 etc b = u 3(x1, x2, ....xN) ... q = u S(x1, x2, ....xN) ( 1 . 1 0 ) then the matrix shown in (1.2) is exactly the R matrix for this general transformation u = F(x). If S < N, then the R matrix is not square, the relation du = Rd x is a "projection" of a larger space into a smaller one, and the equation therefore cannot be inverted. This notation "R" is used extensively in our Tensor Analysis document, see in particular Chapter 2, though only invertible transformations F:EN→ EN are considered there. Sometimes R is referred to as "the differential" of a transformation. Before proving Theorem 1, we want to cut to th e chase and state the Method of Lagrange Multipliers to create a framework in which candidate soluti ons for the extremum problem can be found. 6 2. The Method of Lagrange Multipliers The extrem um problem stated above in (1.1), find extremum of f( r) subject to constraints a( r)=0, b( r)=0, c(r )=0, .... q(r )=0 . (1.1) 1 2 3 S-1 can be solved in the following manner: (a) Define the following function H [ the λ i are the "Lagrange multipliers" ] H( r,λ) ≡ f(r) + λ 1a(r) + λ2 b(r) + ...... + λS-1 q(r) . ( 2 . 1 ) Notice that H is a function of N + (S-1) = N+S-1 variables which are the x i and λi. (b) Compute all N+S-1 partial derivatives of H, which as usual we denote by H i. For the first N partials i = 1,2..N one gets H i = fi + λ1ai + λ2 bi+ ...... + λS-1 qi i = 1,2....N . (2.2) For the final S-1 partials one finds, H N+1 = a HN+2 = b ... H N+S-1 = q . ( 2 . 3 ) (c) Now set all N+S-1 partials of H to 0, that is, set H i = 0 for i = 1 to N+S-1. Then (2.2) and (2.3) read f i(r) + λ1ai(r) + λ2 bi(r) + ...... + λS-1 qi(r) = 0 i = 1,2....N // = H i(r) a( r) = 0 b(r ) = 0 .... q( r) = 0 . i = N+1...N+S-1 (2.4) The second group of S-1 equations is a statement of the S-1 constraints of the problem. The first set of equations in (2.4) can be written in this vector notation f(r) + λ 1a(r) + λ2 b(r) + ...... + λS-1 q(r) = 0 // = H(r) (2.5) where each bolded vector is a row of the R matrix (1.2 ). This equation has the fo rm (1.6) and certainly the coefficient of f(r) is non-zero, being one. Equation (2.5) thus says that the S rows of the R matr ix are linearly dependent! According to (1.7), this means that rank[R(r)] < S. If one had rank(R) = S, all S rows of R would be linear ly independent and then an equation of the form (2.5) could not exist. According to Theorem 1 of (1.8), any point r which makes 7 (2.5) be true is therefore a candidate for a solution to the constrained extremum problem (1.1). Such a solution value of r will be associated with solution values for the λi. (d) It is then incumbent on the problem solver to solve the following N+S-1 equations for the N+S-1 unknowns which are r = ( x1, x2 .....xN) and λ1, λ2.....λS-1 . Here again are the equations to solve: f 1(r) + λ1a1(r) + λ2 b1(r) + ...... + λS-1 q1(r) = 0 f2(r) + λ1a2(r) + λ2 b2(r) + ...... + λS-1 q2(r) = 0 // recall: b 2 means ∂b/∂x2 .... f N(r) + λ1aN(r) + λ2 bN(r) + ...... + λS-1 qN(r) = 0 // N equations a( r) = 0 b( r) = 0 ... q( r) = 0 // S-1 constraint equations (2.6) Since all the functions in these equations are known (since f( r) and all the constraints are known), it is in principle possible to solve the equations for r and the λ i. Then those solution r values are candidate solutions to the constrained extremum problem and need to be studied more to see if they really solve the problem of interest (max, min, saddle, et c). In general the various functions of r are non-linear and not just polynomials, and there are likely to be multip le candidate solutions. A numeric solution may be required. At least this method very significantly re duces the search for solutions to the initial problem (1.1). The Method of Lagrange Multipliers does not provide a silver bullet for solving these equations. Most examples one sees in texts and on the web have very simple functions allowing a straightforward analytic solution. Since the λ i appear linearly in the first N equations in (2 .6), it is possible to obtain expressions for the λi as follows. First, write the first N equations of (2.6) in matrix form, - ⎝⎜⎛ ⎠⎟⎞f1 f2 ... fN = ⎝⎜⎛ ⎠⎟⎞a1 b1 ... q1 a2 b2 ... q2 ... ... ... ... aN bN ... qN ⎝⎜⎛ ⎠⎟⎞λ1 λ2 ... λS-1 ( 2 . 7 ) where the matrix has S-1 columns (one column for each constraint) and N rows. Since S ≤ N, one has S-1 < N so there are fewer columns than rows. This matrix is the transpose of the R matrix (1.2) with the f i elements removed. Taking only the first S-1 rows of the above matrix equation, one finds - ⎝⎜⎛ ⎠⎟⎞f1 f2 ... fS-1 = ⎝⎜⎛ ⎠⎟⎞a1 b1 ... q1 a2 b2 ... q2 ... ... ... ... aS-1 bS-1 ... qS-1 ⎝⎜⎛ ⎠⎟⎞λ1 λ2 ... λS-1 ( 2 . 8 ) 8 where the matrix is now square. Assuming this matrix has non-zero determinant, one can invert (2.8) to obtain expressions for the Lagrange multipliers λi(r) as functions or r. These λ i functions can be inserted into the first N equations of (2.6) and then one can proceed to solve those N equations for r . By making different versions of (2.8) taking diffe rent sets of S-1 rows from (2.7), one may produce several different sets of { λi(r)} corresponding to different candidate solutions r. Stationary Points vs. Extremum Points The Method of Lagrange Multipliers really identifi es stationary points of f some of which may be extremum points. A stationary point (critical point ) of f(x 1,x2, x3....xN) is one where all the partial derivatives vanish. In the notation above, we would write f i = 0 for i = 1,2...N or just f = 0. The Lagrange Multiplier method looks for stationary poi nts of the auxiliary function H where H = 0 as in (2.5). For f(x,y) the reader well knows that a stationary point where ∂xf=0 and ∂yf = 0 could be a saddle point which is neither a minimum nor a maximum. There could be many stationary points for a given problem, some of which are extrema, and some of which are not ( see Example 3 in Section 5). Each stationary point will have its own set of Lagrange Multipliers { λi}. Moreover, a true extremal point may occur at a boundary of the domain of interest, and would then not be detected by the Lagrange Multiplier method. The following 1D function f(x) illustrates many of these points. ( 2 . 9 ) Here the domain of interest is the portion of the x axis lying between the two vertical bars. Point a is the true minimum and occurs at a boundary. Point b is the true global maximum, while point e is only a local maximum. Point d wants to be the global minimum but is trumped by point a. Point c has ∂ xf = 0 but is neither a local minimum nor a local maximum -- it is the 1D analog of a saddle point. Had points b and e the same height, there would be two equal global maxi ma. All these situations can occur for the function H for the constrained extremum problem and so analysis is required after the stationary points of H are located. Comments 1. For the Lagrange Multiplier method to be valid, the functions f( r), a(r), b(r)....q( r) must be continuous and differentiable in all arguments x i and their derivatives must also be continuous. The functions are therefore C0 and C1. This ensures that the R( r) matrix in (1.2) behaves smoothly (is continuous) as a function of r. 9 2. Notation : Some authors use the notation f xi≡ ∂f/∂xi and the constraints are g(s)(x) = 0 where s is a label, Then the first derivative of a constraint function is g(s) xi(x). Another choice is to refer to the constraint functions as g s(x) where s is a label, and then one just writes out ∂gs ∂xi . We feel that our simple though less rigorous notation provides more clarity. Writing ∂f/∂xi = fi is a notation promoted by Buck. 3. The λi terms in (2.1) often appear with minus signs in place of plus signs. This removes the minus sign in (2.7) and (2.8). These signs are just a convention and have no effect on the solution points r. 4. The function H in (2.1) is sometimes referred to as the Lagrangian and is written as L. This Lagrangian differs from that which appears in the Lagrangian formulation of particle motion. In Lagrangian dynamics, the motion of a set of particles (perhaps comprising a rigid body) can be determined from a certain Lagrangian L ≡ T - V where T and V are the kinetic and potential energies of the set of particles. L is a function of the "generalized coordinates" q k and velocities q •k = dqk/dt of the particles, and of time t, so one writes L(q k,q•k, t). The equations of motion in the absence of constraints are known as the Euler-Lagrange equations, ( ∂L ∂qk – d dt ∂L ∂q•k ) = 0. One may add a set of differential constraints to the problem of the form Σkaikdqk = 0 where the a ik are functions causing the particles to be constrained in their motion. For example, the pa rticles may be constrained to move only along certain curves like beads on a wire. The adjusted Euler-Lagrange equations then have the form ( ∂L ∂qk – d dt ∂L ∂q•k ) + Σiλiaik = 0. This equation is similar to the first line of our equation (2.4). See Goldstein page 41 (2.23) and (2.27). Goldstein refers to this method of incor porating constraints into particle dynamics as the method of Lagrange undetermined multipliers. 10 3. Proof of Theorem 1 Theore m 1. If the point r is a solution to the extremal problem (1.1), then rank[R( r)] < S. (1.8) (a) Proof for N = 6 and S = 4 The more general the proof is made, the less clear it becomes due to the notational baggage that must be added. Therefore we shall prove Theorem 1 for the specific case shown above in Fig (1.3) where N = 6 and S = 4 (3 constraints). The reader should have no trouble following each step for this sample case. When we reach the proof conclusion for the sample case, we start over again doing the general case and compare each step to the steps of this sample case. Although a lot of e quations appear below (and a lot of column inches are consumed), everything is really quite simple. So, for our sample case we start off with this R ma trix, again with bolded row vectors on the right, f 1 f2 f3 f4 f5 f6 f a 1 a2 a3 a4 a5 a6 a R = b 1 b2 b3 b4 b5 b6 = b (3.1) c 1 c2 c3 c4 c5 c6 c We want then to find an extremum of f(x 1, x2, x3, x4, x5, x6) subject to these constraints: a(x 1, x2, x3, x4, x5, x6) = 0 b(x1, x2, x3, x4, x5, x6) = 0 c(x1, x2, x3, x4, x5, x6) = 0 . ( 3 . 2 ) In practice, a person solving this problem might try to eliminate variables. For example, somehow with enough labor one should be able to "use up" the three c onstraint equations to eliminate three of the six variables x i. There are of course (6,3) ways to select three variables for elimination. For starters, imagine that we have used up the 3 constraints to eliminate variables x 1, x2 and x3. Then we write a(x 1, x2, x3, x4, x5, x6) = 0 x 1 = X1(x4, x5, x6) b(x1, x2, x3, x4, x5, x6) = 0 ⇒ x2 = X2(x4, x5, x6) c(x1, x2, x3, x4, x5, x6) = 0 x 3 = X3(x4, x5, x6) (3.3) where X1, X2 and X3 are three resulting functions. This is a "theoretical" elimination, one does not actually have to do the process, one need only underst and that in principle it could be done and the functions Xi could in principle be found. Now rewrite the function f and the three constraint f unctions in this manner, substituting for example the function X1 for x1 everywhere it appears, 11 f(X1(x4, x5, x6), X2(x4, x5, x6), X3(x4, x5, x6), x4, x5, x6) ≡ F(x4, x5, x6) a(X1(x4, x5, x6), X2(x4, x5, x6), X3(x4, x5, x6), x4, x5, x6) ≡ A(x4, x5, x6) b(X1(x4, x5, x6), X2(x4, x5, x6), X3(x4, x5, x6), x4, x5, x6) ≡ B(x4, x5, x6) c(X1(x4, x5, x6), X2(x4, x5, x6), X3(x4, x5, x6), x4, x5, x6) ≡ C(x4, x5, x6) . (3.4) In this way we have defined four new functions F,A,B,C each of the 3 variables shown (those that were not eliminated). Next, compute the first partial derivatives of F using the chain rule. For example ∂F ∂x4 = ∂f ∂x1 ∂X1 ∂x4 + ∂f ∂x2 ∂X2 ∂x4 + ∂f ∂x3 ∂X3 ∂x4 + ∂f ∂x4 . ( 3 . 5 a ) In our compact notation where F i means ∂F/∂xi this can be restated as F4 = f1X1 4 + f2X2 4 + f3X3 4 + f4 . ( 3 . 5 b ) Doing this for each of the non-eliminated variables one gets, F4 = f1X1 4 + f2X2 4 + f3X3 4 + f4 F5 = f1X1 5 + f2X2 5 + f3X3 5 + f5 F6 = f1X1 6 + f2X2 6 + f3X3 6 + f6 . ( 3 . 6 ) Notice that the subscripts on the leftmost three f i correspond to those of the eliminated variables i = 1,2,3. Now, since we are seeking an extremal point of F(x 4, x5, x6) with no constraints at all (the three initial constraints have all been absorbed in the process of eliminating three variables), we require that (x 4, x5, x6) be a "critical point" for the functio n F, which means we require that F 4 = F5 = F6 = 0. Then the above equations become, f1X1 4 + f2X2 4 + f3X3 4 + f4 = 0 f1X1 5 + f2X2 5 + f3X3 5 + f5 = 0 f1X1 6 + f2X2 6 + f3X3 6 + f6 = 0 . ( 3 . 7 ) Now apply the same process to the three constraint functions A,B,C. First compute the derivatives, and then set the derivatives to zero. One does this for A(x 4, x5, x6), for example, because the constraint condition a = 0 says that that A(x 4, x5, x6) = 0 = a constant, for all values of x 4, x5, x6 , so all first derivatives must vanish. Doing this for f unction A then gives these three equations, a1X1 4 + a2X2 4 + a3X3 4 + a4 = 0 a1X1 5 + a2X2 5 + a3X3 5 + a5 = 0 a1X1 6 + a2X2 6 + a3X3 6 + a6 = 0 . ( 3 . 8 ) We have just taken the previous equations and replaced f → a everywhere. Then do this for B and C. One ends up then with this set of 12 equations: 12 f1X1 4 + f2X2 4 + f3X3 4 + f4 = 0 f1X1 5 + f2X2 5 + f3X3 5 + f5 = 0 f1X1 6 + f2X2 6 + f3X3 6 + f6 = 0 a1X1 4 + a2X2 4 + a3X3 4 + a4 = 0 a1X1 5 + a2X2 5 + a3X3 5 + a5 = 0 a1X1 6 + a2X2 6 + a3X3 6 + a6 = 0 b1X1 4 + b2X2 4 + b3X3 4 + b4 = 0 b1X1 5 + b2X2 5 + b3X3 5 + b5 = 0 b1X1 6 + b2X2 6 + b3X3 6 + b6 = 0 c1X1 4 + c2X2 4 + c3X3 4 + c4 = 0 c1X1 5 + c2X2 5 + c3X3 5 + c5 = 0 c1X1 6 + c2X2 6 + c3X3 6 + c6 = 0 . ( 3 . 9 ) Here there are S = 4 equation groups, and each group h as N-S+1 = 6-4+1 = 3 equations (this last 3 is the number of non-eliminated variables). Next reorder these equations into three groups of four, for example taking the first equation in each group above to make the first new group below, f 1X1 4 + f2X2 4 + f3X3 4 + f4 = 0 a1X1 4 + a2X2 4 + a3X3 4 + a4 = 0 b1X1 4 + b2X2 4 + b3X3 4 + b4 = 0 c1X1 4 + c2X2 4 + c3X3 4 + c4 = 0 f 1X1 5 + f2X2 5 + f3X3 5 + f5 = 0 a1X1 5 + a2X2 5 + a3X3 5 + a5 = 0 b1X1 5 + b2X2 5 + b3X3 5 + b5 = 0 c1X1 5 + c2X2 5 + c3X3 5 + c5 = 0 f 1X1 6 + f2X2 6 + f3X3 6 + f6 = 0 a1X1 6 + a2X2 6 + a3X3 6 + a6 = 0 b1X1 6 + b2X2 6 + b3X3 6 + b6 = 0 c1X1 6 + c2X2 6 + c3X3 6 + c6 = 0 . (3.10) Now there are N-S+1 = 3 equation groups , and each group has S = 4 equations. Next, write each of these three equation sets as a matrix equation, 13 ⎝⎜⎛ ⎠⎟⎞ 0 0 0 0 = ⎝⎜⎛ ⎠⎟⎞f1 f2 f3 f4 a1 a2 a3 a4 b1 b2 b3 b4 c1 c2 c3 c4 ⎝⎜⎛ ⎠⎟⎞ X1 4 X2 4 X3 4 1 = (f 1,f2, f3, f4) ⎝⎜⎛ ⎠⎟⎞ X1 4 X2 4 X3 4 1 ⎝⎜⎛ ⎠⎟⎞ 0 0 0 0 = ⎝⎜⎛ ⎠⎟⎞f1 f2 f3 f5 a1 a2 a3 a5 b1 b2 b3 b5 c1 c2 c3 c5 ⎝⎜⎛ ⎠⎟⎞ X1 5 X2 5 X3 5 1 = (f 1,f2, f3, f5) ⎝⎜⎛ ⎠⎟⎞ X1 5 X2 5 X3 5 1 ⎝⎜⎛ ⎠⎟⎞ 0 0 0 0 = ⎝⎜⎛ ⎠⎟⎞f1 f2 f3 f6 a1 a2 a3 a6 b1 b2 b3 b6 c1 c2 c3 c6 ⎝⎜⎛ ⎠⎟⎞ X1 6 X2 6 X3 6 1 = (f 1,f2, f3, f6) ⎝⎜⎛ ⎠⎟⎞ X1 6 X2 6 X3 6 1 . (3.11) On the right we show our abbreviated notation for a matrix, just showing the first row. Now we claim that the three matr ices shown must have zero determinant! Suppose the first matrix equation had a non-zero determinant. It could then be inverted (see (B.4.8)) to give ⎝⎜⎛ ⎠⎟⎞ X1 4 X2 4 X3 4 1 = (f 1,f2, f3, f4)-1 ⎝⎜⎛ ⎠⎟⎞ 0 0 0 0 = ⎝⎜⎛ ⎠⎟⎞ 0 0 0 0 . (3.12) But this produces a blatant contradiction that 1 = 0 (not to mention the other rows), and therefore we must have det (f 1,f2, f3, f4) = 0 det (f 1,f2, f3, f5) = 0 det (f 1,f2, f3, f6) = 0 or det (f 1,f2, f3, fm) = 0 m ≠ 1,2,3 . (3.13) Therefore, we have shown that three of the SxS = 4x4 submatrices of Fig (1.3) vanish. But there are 15 such submatrices, and in order to show th at rank(R) < S, we have to show that all 15 SxS determinants vanish. But that can be demonstrated by starting over several times and each time c hoosing to eliminate three different variables. Notice in the det = 0 equations above that the first three f i values correspond to the three eliminated variables x 1, x2 and x3, Had we instead chosen to eliminate x 1, x2 and x4, we would have obtained, 14 det (f 1,f2, f4, f3) = 0 det (f 1,f2, f4, f5) = 0 det (f 1,f2, f4, f6) = 0 or det (f 1,f2, f4, fm) = 0 m ≠ 1,2,4 . (3.14) More generally, had we chosen to eliminate variables x i, xj and xk ( i ≠ j ≠ k), we would have obtained det (f i,fj, fk, fm) = 0 m ≠ i , j , k . ( 3 . 1 5 ) As i,j,k range over all possible values 1,2,3,4,5,6, we clearly hit all possible 4x4 subdeterminants, and thus we have shown that they all vanish, and therefore rank(R) < S, QED. Lest it be overlooked, one must keep in mind that r = (x1, x2, x3, x4, x5, x6) must be a solution to the extremum problem (or at least a critical point as in (3 .7) ) including meeting all the constraints. And then the conclusion is that, for such a solution point r, one has rank[R( r)] < S . Having now shown that rank[R( r)] < S for a solution point r, we know from (1.7) (b) that the S rows of the R matrix must be linearly dependent for that r . This means that one can find a set of constants λ1,λ2,λ3 such that f(r) + λ1a(r) + λ2 b(r) + λ3 c(r) = 0 . (2.5) These constants are the Lagrange Multipliers and we see that the fact that they exist follows directly from the fact that rank[R( r)] < S. To summarize, one theoretically eliminates N-S+1 of the N variables to create the functions F,A,B,C of (3.4). One then looks for the normal (unconstrained) critical points of F where all partials vanish. Since the constraints all have the form A = 0, B = 0, C = 0, their partials vanish as well. Writing all these equations in matrix form (for various sets of eliminat ed variables) leads to the conclusion that all the SxS minors of matrix R vanish so that rank(R) < S. (b) Proof for general S ≤ N To generalize now t o general N and S, we restate some of the previous equations (adding a prime to the equation number) using a more generalized notation. The reader w ill note how even our simple notation becomes rather unpleasant. Comment : This complexity increases further if one denotes the constraint functions by a notation like a = g(1), b = g(2) .... q = g(S-1). We avoid the idea of a = g 1, b = g2 etc because we want a subscript to indicate a partial derivative, not a label. We start off with this R matrix which has N columns and S ≤ N rows. The S-1 constraint functions are a,b,....q. Remember that f 2 means ∂ f/∂x2. 15 f 1 f2 f3 ... f N f a 1 a2 a3 ... a N a R = b 1 b2 b3 ... b N = b (3.1)' .... ... q 1 q2 q3 ... q N q Here are the S-1 constraint equations, a(x 1, x2, ...xN) = 0 // S-1 constraint equations b(x1, x2, ...xN) = 0 ... q(x 1, x2, ...xN) = 0 . ( 3 . 2 ) ' Eliminate (theoretically) the first S-1 variables x 1,x2....xS-1 using the S-1 constraint equations: a(x1, x2, ...xN) = 0 x 1 = X1(xS, xS+1 ... xN) b(x1, x2, ...xN) = 0 ⇒ x2 = X2(xS, xS+1 ... xN) ... ... q(x 1, x2, ...xN) = 0 x S-1 = XS-1(xS, xS+1 ... xN) . (3.3)' Substitute to get the F,A,B...Q equations (a set of S equations), f(X 1(xS, xS+1 ... xN), X2(xS, xS+1 ... xN), .. XS-1(xS, xS+1 ... xN) , xS, xS+1 ... xN) ≡ F(xS, xS+1 ... xN) a(X1(xS, xS+1 ... xN), X2(xS, xS+1 ... xN), .. XS-1(xS, xS+1 ... xN) , xS, xS+1 ... xN) ≡ A(xS, xS+1 ... xN) b(X1(xS, xS+1 ... xN), X2(xS, xS+1 ... xN), .. XS-1(xS, xS+1 ... xN) , xS, xS+1 ... xN) ≡ B(xS, xS+1 ... xN) ... q(X 1(xS, xS+1 ... xN), X2(xS, xS+1 ... xN), .. XS-1(xS, xS+1 ... xN) , xS, xS+1 ... xN) ≡ Q(xS, xS+1 ... xN) . ( 3 . 4 ) ' Now take first derivatives as shown in (3.5) to get this generalized version of (3.6), F S = f1X1 S + f2X2 S + f3X3 S + ... + f S-1XS-1 S + fS FS+1 = f1X1 S+1 + f2X2 S+1 + f3X3 S+1 + ... + f S-1XS-1 S+1 + fS+1 FS+2 = f1X1 S+2 + f2X2 S+2 + f3X3 S+2 + ... + f S-1XS-1 S+2 + fS+2 ... F N = f1X1 N + f2X2 N + f3X3 N + ... + f S-1XS-1 N + fN . // N-S+1 equations (3.6)' Requiring (x S, xS+1 ... xN) to be a "critical point", we set all the these first derivatives to 0 to get, 16 f1X1 S + f2X2 S + f3X3 S + ... + f S-1XS-1 S + fS = 0 f1X1 S+1 + f2X2 S+1 + f3X3 S+1 + ... + f S-1XS-1 S+1 + fS+1 = 0 f1X1 S+2 + f2X2 S+2 + f3X3 S+2 + ... + f S-1XS-1 S+2 + fS+2 = 0 ... f1X1 N + f2X2 N + f3X3 N + ... + f S-1XS-1 N + fN = 0 . // N-S+1 equations (3.7)' A similar set of equations is obtained for a,b,c...q. Just replace f →a, then f → b and so on. We shall not write out all these equations as we did in (3.9). We have then S sets of equations, each set containing N-S+1 equations. As in the sample case, we then re order the equations to get first this group, f 1X1 S + f2X2 S + f3X3 S + ... + f S-1XS-1 S + fS = 0 a1X1 S + a2X2 S + a3X3 S + ... + a S-1XS-1 S + aS = 0 b1X1 S + b2X2 S + b3X3 S + ... + b S-1XS-1 S + bS = 0 .... q 1X1 S + q2X2 S + q3X3 S + ... + q S-1XS-1 S + qS = 0 . // S equations (3.10)' S The next group has S →S+1 in all subscript positions. The last group has S →N, and here is that last group, f 1X1 N + f2X2 N + f3X3 N + ... + f S-1XS-1 N + fN = 0 a1X1 N + a2X2 N + a3X3 N + ... + a S-1XS-1 N + aN = 0 b1X1 N + b2X2 N + b3X3 N + ... + b S-1XS-1 N + bN = 0 .... q 1X1 N + q2X2 N + q3X3 N + ... + q S-1XS-1 N + qN = 0 . // S equations (3.10)' N So there are now N-S+1 groups of equations and each group has S equations each having S terms. The next task is to write each set of S equations as a matrix equation. Here we do that just using the abbreviated notation for the matrix, ⎝⎜⎜⎛ ⎠⎟⎟⎞ 0 0 ... 0 0 = (f 1,f2, f3, ... fS-1, fS) ⎝⎜⎜⎛ ⎠⎟⎟⎞ X1 S X2 S ... XS-1 S 1 ⎝⎜⎜⎛ ⎠⎟⎟⎞ 0 0 ... 0 0 = (f 1,f2, f3, ... fS-1, fS+1) ⎝⎜⎜⎛ ⎠⎟⎟⎞ X1 S+1 X2 S+1 ... XS-1 S+1 1 •••• ⎝⎜⎜⎛ ⎠⎟⎟⎞ 0 0 ... 0 0 = (f 1,f2, f3, ... fS-1, fN) ⎝⎜⎜⎛ ⎠⎟⎟⎞ X1 N X2 N ... XS-1 N 1 . ( 3 . 1 1 ) ' 17 We then argue as before that, to avoid the contradic tion 0 = 1, the determinants of these SxS submatrices must be zero, det(f 1,f2, f3, ... fS-1, fS) = 0 det(f 1,f2, f3, ... fS-1, fS+1) = 0 ... det(f 1,f2, f3, ... fS-1, fN) = 0 or det(f 1,f2, f3, ... fS-1, fm) = 0 m ≠ 1,2,3...S-1 . (3.13)' Finally, if we start off by eliminating the S-1 variables x i, xj, xk.... xr where i ≠ j ≠k....≠ r, we would find that det(f i,fj, fk, ... fr, fm) = 0 m ≠ i,j,k...r . (3.15)' Since this includes any possible SxS subdeterminant of the R matrix, we have then shown that all (N,S) SxS subdeterminants vanish, and therefore rank[R( r)] < S when r is a solution or critical point of the extremum problem. Having now shown that rank[R( r)] < S for a solution point r, we know that the S rows of the R matrix must be linearly dependent for that r . This means that one can find a set of constants λ 1,λ2,...λS-1 such that f(r) + λ 1a(r) + λ2 b(r) + ...... + λS-1 q(r) = 0 . (2.5) These constants are the Lagrange Multipliers and we see that the fact that they exist follows directly from the fact that rank[R( r)] < S. 18 4. The Gradient Interpretation: Examples 1 and 2 4.1 Setup Recall that, for r such that rank[R( r)] < S, one can find Lagrange multipliers λi such that f(r) + λ 1a(r) + λ2 b(r) + ...... + λS-1 q(r) = 0 , (2.5) where the bolded letters are the rows of the R matrix shown in (1.2). In component notation this equation reads fi(r) + λ1ai(r) + λ2 bi(r) + ...... + λS-1 qi(r) = 0 i = 1,2...N (4.1.1) where recall f i = ∂f/∂xi. In the usual gradient operator ∇ notation this can be written ∇f(r) + λ1∇a(r) + λ2 ∇b(r) + ...... + λS-1 ∇q(r) = 0 ∇ = (∂1,∂2...∂N) . (4.1.2) Having obtained this equation by the matrix rank met hod, we shall now show that it makes perfect sense in terms of the gradients. We shall consider two Exam ples below in order to build up to the general claim made at the end in Section 4.4. The reader unint erested in examples may skip to that section. Before continuing, here are a few Facts which perhaps are "old hat" to the reader but need stating: Fact 1: The equation F(x 1,x2.....xn) = K defines an n-1 dimensional surface in En . (4.1.3) For example, F(x,y,z) = x2+y2+z2 with F = R2 describes a spherical 2D surface of radius R in E3. F(x,y) = x2+y2 with F= R2 describes a 1D surface (curve) in E2, a circle Fact 2: ∇F(r) is always locally normal to the surface F(r ) = K, where r = (x1,x2.....xn) . (4.1.4) Proof : Start at point r on the surface F( r) = K and move a small distance d r in an arbitrary direction along the surface. Since one stays on the surface, F( r+dr) = K. Then dF = F( r+dr)- F( r) = 0. But one knows that dF = ∇F • dr so, in order that dF = 0 for all dr displacements on the surface, ∇F(r) must be locally normal to the surface at point r. Fact 3: As K takes a set of different values K 1, K2.....KS, the equation F(x 1,x2.....xn) = K describes a family of so-called level surfaces (level sets) . If S is large and the range of the K i is small, these surfaces will be closely spaced. Adjacent su rfaces then have nearly the same shape, although in general all the surfaces in the family of surfaces do not have the same shape. (4.1.5) 19 Example : F(x,y) = x2+2y4 and F= K for K = 1 to 10. Here "level surfaces" means "level curves". (4.1.6) Fact 4: The gradient ∇F points in the direction in which F changes most rapidly. (4.1.7) Proof : For tiny displacement d r, dF = ∇F • dr is largest positive when d r points in the direction of ∇ F ("uphill"). And dF = ∇ F • dr is largest negative when d r points opposite the direction of ∇F ("downhill"). With these Facts established, we now consider a simple Lagrange Multiplier example. 4.2 Example 1: Half sphere with one simple constraint Let u = f(x,y) represent the surface of a sphere (radius R = 2) in E3, centered at the origin. We consider only the upper half of this surface, and we want to find (x,y) that maximizes f. If there are no constraints, then (4.1.2) above says ∇f = 0. The only place on our u = f surface having ∇f = 0 is the north pole of the sphere. This is just a regular "critical point" where ∂ xf = 0 and ∂yf = 0, so we are happy with this interpretation of (4.1.2) with no constraints. We now add a constraint a(x,y) = 0 where a(x,y) = y- 1. This constraint is the line y=1 in the x-y plane E2. We can extrude this line into a plane y = 1 in E3. The hemispherical surface is a 2D surface in E3, and the extruded constraint pl ane is also a 2D surface in E3. These 2D surfaces intersect in a 1D surface which is a curve. There is hopefully some point on this curve that maximizes f. Below is a picture of the sphere. The intersection of the upper spherical surface with the plane y = 1 is shown as a red curve (a half circle ). Due to the constraint, we cannot get to the north pole so the maximum value of f is some value less than that u = 2 at the north pole. The extr emum of this constrained problem will be at point A , and point B is not an extremum. 20 (4.2.1) What is potentially misleading in this example (so far) is the notion of surface and gradient. For the above sphere, the surface is described by r = 2 in spherical coordinates, and then if f(r ) = r, one finds that ∇f = ∂rf r^ = r^, and this is indeed normal to the surface at ev ery point on it, as predicted by Fact 2 above. Similarly, the y=1 plane constraint function a(x,y) = y-1 = 0 has ∇a = y^ and this is everywhere normal to the constraint plane. The potential confusion is that the Lagrange Multiplier Show does not play out on the stage shown in (4.2.1). It plays out in a space of one lower dime nsion. For Example 1, the Lagrange Multiplier Show plays out on the disk which lies at u = 0 in (4.2.1). The gradients involved with the Lagrange Multiplier method for Example 1 are 2D gradients which lie in th is disk, not the 3D gradients to surfaces appearing in (4.2.1). So let us move then to the proper setting which is th is u=0 disk in the x,y plane: (viewed from above) (4.2.2) 21 We wish to maximize the function f(x,y) = 4 - x2- y2 subject to the constraint y = 1. In (4.2.1) the spherical surface is u2+x2+y2 = R2 = 4, so f(x,y) is then height u. So the function 4 - x2- y2 = 4 - r2 is the height of the hemisphere in Fig (4.2.1) above th e point (x,y). If we consider f(x,y) = K for a set of K values, we get a set of level curves in (4.2.2) which are circles with r = 4-K . On each circle in (4.2.2), the height of the sphere lying over the surface in (4.2 .1) is constant -- that is why they are called level curves. The constraint is a(x,y) = 0 with constraint functi on a(x,y) = y-1, and the constraint y = 1 is shown as the red line in (4.2.2). We compute, f = 4 - x2- y2 = 4 - r2 a = y - 1 ( 4 . 2 . 3 ) so ∇f = ∇4 - r2 = [-r/ 4-r2 ] r^ ∇a = ∇ ( y-1) = y^ . ( 4 . 2 . 4 ) For any point r in the disk of (4.2.2), ∇f therefore points toward the disk center, and this is then the uphill direction for the hemisphere in (4.2.1). In Cartesian coordinates, ∇f = ∇4 - x2- y2 = - (x/f) x^ - (y/f) y^ = - x 4 - x2- y2 x^ - y 4 - x2- y2 y^ . (4.2.5) Now, equation (4.1.2) says ∇f(r) = – λ1∇a(r) = 0 . ( 4 . 2 . 6 ) How do we make geometric sense of this equation? If r is an "extremum" location, we should have df = 0 for any legal small displacement d r starting at this point r. Eq. (4.2.6) dotted into d r then says, 0 = df = ∇ f(r) • dr = – λ1∇a(r) • dr . (4.2.7) If we displace d r in any direction on the constraint surface (which from Fact 2 has normal ∇a(r)), then ∇a(r) • dr = 0 and we find df = 0, so we are thus at an extremum of f. For other values of position r , we will find either that df > 0 or df < 0 if we move d r along the constraint surface. If we are looking for a maximum of f, for example, then it is to our advantage to move a little d r in a direction for which df > 0. Doing this repeatedly to find a maximum (minimum) of f is called the method of gradient ascent (descent ) In our example with only one constraint, (4.2.6) says that we have reached an extremum when ∇f and ∇a are collinear. Let's draw onto (4.2.2) the points a and b which lie under points A and B in (4.2.1), and we include a point b' which is a mirror image of b : 22 (4.2.8) At each of the points a,b,b' we show ∇f in black and ∇a in red. At point a we have reached the closest distance to the center that is allowed for points on th e red constraint line, so this will correspond to the maximum value of f, which is point A in (4.2.1). At point a we see that in fact the two gradients are collinear as required by (4.1.2) at an extremum. Suppose we are at point b. We may compute df for a small displ acement upwards (minus x direction) df( b) = ∇ f • dr = [ - (x/f) x^ - (y/f) y^ ] • |dx|(- x^) = |dx| (x/f) > 0 since x>0 and f>0 . (4.2.9) Since df > 0, it is to our advantage in finding max f to move upwards in (4.2.8). Conversely, suppose we are instead at point b'. Then moving toward a gives, df( b') = ∇ f • dr = [ - (x/f) x^ - (y/f) y^ ] • |dx|( x^) = |dx| (-x/f) > 0 since x<0 and f>0 . (4.2.10) and again we find that df > 0. From either starting position b or b', moving toward point a is a win. It is an easy matter to compute the so lution value of the Lagrange Multiplier λ 1 and r. Inserting (4.2.4) into (4.2.6) gives ∇f(r) = – λ 1∇a(r) ⇒ - x 4 - x2- y2 x^ - y 4 - x2- y2 y^ = -λ1 y^ . (4.2.11) This says x = 0, and then using the constraint y = 1, - 1 4 - 1 y^ = -λ1 y^ ⇒ λ1 = 1/ 3 . (4.2.12) The solution to the Example 1 extremum problem is then 23 r = (0,1) λ1 = 1/ 3 . ( 4 . 2 . 1 3 ) We have claimed in (2.1) that at the solution poin t, the Lagrangian function H(x,y) should have a maximum, since the whole theory is based on H being the function to maximize without constraints. One has, H(x,y) = f(r ) + λ 1a(r) = 4 - x2- y2 + (1/ 3 ) (y-1). (4.2.14) Now u = (1/ 3 ) (y-1) is a plane sloping up to the right in Fig (4.2.1). This is different from the plane y = 1 which is the vertical constraint plane in that figur e. In (4.2.14) we are adding a spherical surface to a plane sloping up to the right, and the result is an e llipsoidal-like surface (really a quartic surface) which has a maximum at the point A = (0,1). Here is a plot of that surface: (4.2.15) One can see from these plots that the maximum of H(x,y) occurs at x = 0 and y = 1. A direct method of confirming the maximum is provided by examining H ii : Hi = fi + λ1ai H ii = fii + λ1aii a ii = ∂2(y-1)/∂ xi2 = 0 f = 4 - x2- y2 = 4 - x12- x22 ⇒ fi = - xi/f f > 0 fii = - [f * 1 - x ifi] / f2 = - [f - x i(-xi/f)] / f2 = - (1/f) - (x i/f)2 H ii = fii + λ1aii = fii = - (1/f) - (x i/f)2 < 0 . (4.2.16) 24 This shows that the function H(x,y) is "cupping down " at all locations, as the graph suggests. The point A where H i = 0 (see (2.4) ) thus also has H ii < 0 and is thus a maximum. If we had a surface in (4.2.1) more complicated than a sphere, and a constraint more complicated than y = 1, the nature of the interpretation of (4.1.2) does not change. For example, if (4.2.1) contained a surface u = f(x,y) whose level curves were those of (4.1.6) with f = K values increasing toward the center, and if the constraint a(x,y) = 0 were some arbitrary constraint curve (shown in red), we would have this picture: (4.2.17) In this case the red arrow ∇a changes direction as one moves along the red curve, but the solution point is shown as the black dot a where the red arrow and the black arrow representing ∇f are collinear. At point b the gradient arrows do not line up, and there is advantage for increasing f by moving toward point a. Notice that saying the two arrows are collinear is to sa y they are linearly dependent. Equation (4.1.2) says that the gradients are linearly dependent in the general case as well as this simple case. 4.3 Example 2: Half 4D sphere with two simple constraints In order to better demonstrate the interpretation of (4.1.2) concerning the gradients. we upgrade Example 1 and then allow two constraints. We take the surface of (4.2.1) to be the "upper half" of a 4D sphere, whatever that means. Since we cannot draw such a thing in E4 even in a projection drawing like (4.2.1), we don't even try. However, the Lagrange Multipli er scenario (the analog of (4.2.2)) exists in E3, not E4, and we can at least draw that. In Example 1 the le vel curves were a set of concentric circles in E2. For Example 2 the level surfaces are a set of concentric spheres in E3 each labeled by a value of K, where f(x,y,z) = K and where f(x,y,z) = 4 - x2- y2- z2 . We are not proving this claim, it is just based on analogy: if the slice of a 3D sphere is a 2D sphere (a disk), then a s lice of a 4D sphere ought to be a 3D sphere. We take as our two constraints the equations y = 1 and x = 1, just to keep it simple. Thus 25 f(x,y,z) = 4 - x2- y2- z2 = 4 - r2 a(x,y,z) = y-1 b ( x , y , z ) = x - 1 . ( 4 . 3 . 1 ) The gradients of interest are now 3D gradients in E 3 space, so ∇f = ∇ (4 - r2 ) = [-r/ 4-r2 ] r^ = - (x/f) x^ - (y/f) y^ - (z/f) z^ points to sphere center ∇a = ∇ ( y-1) = y^ ∇b = ∇ ( x-1) = x^ . ( 4 . 3 . 2 ) Equation (4.1.2) reads ∇f(r) + λ1∇a(r) + λ2 ∇b(r) = 0 . ( 4 . 3 . 3 ) This says that the three gradients must be linearly de pendent, which means they must be coplanar! Before solving the problem, we would like to draw a picture corresponding to (4.2.2) for Example 1. Although one could in fact draw such a picture with some effort , we shall decline this task and instead draw two z = constant slices of the desired picture. On the left belo w is a slice at z = 0, while on the right is a slice at z = 1 (these are not precision pictures) : ( 4 . 3 . 4 ) In these pictures, we are viewing the two red constraint planes y = 1 and x = 1 edge on. The surface which satisfies both constraints is a line normal to the plane of paper which is th e intersection of the two constraint planes. The circles on the left are slices of the family of concentric spheres taken at the equator z = 0. At z = +1 since our slicing plane is closer to the north pole, some of the inner spheres no longer cut the plane z = 1, and those that still do have smaller diameter circles. 26 Consider now the point r = b shown as a black dot on the right. The two red arrows are ∇a = y^ and ∇b= x^, each normal to its constraint plane. The black arrow dips into the plane of paper because ∇f = [-r/ 4-r2 ] r^ which points toward the sphere center. The three arrows ∇ f, ∇a,∇b are therefore not coplanar, so they are linearly independent. That mean s that equation (4.3.3) cannot exist [ see (1.6) ] for this point r. This black point b is located outside the r = 7 units sphere (in our crude picture). Consider next the point r = a shown as a black dot on the left. Since this is the equatorial slice, the black arrow ∇f lies in the plane of paper so the three arrows are coplanar. This means the three vectors ∇f, ∇a,∇b are coplanar so they are linearly dependent . For this point a, the equation (4.3.3) can and does exist, and therefore point a is the problem solution. This point is located inside the r = 6 unit sphere. Once again, at an extremum point we should have df = 0. Moving an amount d r which is consistent with both constraints (meaning d r is in the z direction) we then have from (4.3.3), df = ∇f • dr = – λ1[∇a(r)• dr] – λ2 [∇b(r)• dr] = –λ1 [0] - λ2[0] = 0 (4.3.5) and sure enough, df = 0. At point b one finds for a d r pointed down toward the equatorial plane, df = ∇f • dr = [ - (x/f) x^ - (y/f) y^ - (z/f) z^ ] • (-|dr| z^ ) = (z/f)dr > 0 (4.3.6) and so it is advantageous to move d r = |dr| z^ toward a and thereby increase f, so b is not an extremum. Evaluating (4.3.6) at z = 0 for point a again shows df = 0. Finally, we solve the problem. Inserting (4.3.2) into (4.3.3) gives ∇f(r) = – λ1∇a(r) – λ2 ∇b(r) ⇒ [- x 4 - x2- y2- z2 x^ - y 4 - x2- y2- z2 y^ - z 4 - x2- y2- z2 z^ ] = – λ1y^ – λ2 x^ . (4.3.7) We see at once that z = 0 and then the above becomes x 4 - x2- y2 = λ2 y 4 - x2- y2 = λ1 . (4.3.8) But the constraints say x = 1 and y = 1 so, 1 4 - 12- 12 = λ2 1 4 - 12- 12 = λ1 ⇒ λ1 = λ2 = 1/ 2 . (4.3.9) Therefore the solution to Example 2 is this: r = (1,1,0) λ 1 = 1/ 2 λ2 = 1/ 2 . (4.3.10) It is simple matter conceptually to generalize our Example 2 to more complicated functions f,a,b. For some arbitrary function f(r ) there will be a set of 2D surfaces in E3 which are the level surfaces on which f(r) = K for various values of K. These level surfaces th en replace the set of concentric spheres of our simple case. We then imagine replacing our simple plan ar constraint functions a(x,y,z) and b(x,y,z) with 27 general functions which then result in arbitrarily shap ed 2D constraint surfaces a(x,y,z) = 0 and b(x,y,z) = 0 in E3. The gradients ∇a and ∇ b vary over their respective surfaces (being normal vectors). The intersection of these two 2D constraint surfaces in E3 will be a curve in E3 (1D surface in E3). Potential solution points r must lie on this curve. At each point r on this curve, the gradients ∇a(r) and ∇b(r) define a local plane. At a stationary point one will find that the gradient ∇f(r) lies in the same plane as ∇a(r) and ∇b(r) so the three gradients are linearly dependent and then ∇f(r) + λ1∇a(r) + λ2 ∇b(r) = 0 for some λ1 and λ2. At such a point r, if one considers any d r which lies on both constraint surfaces, one will then find that df = ∇f •dr = 0 so the point r is therefore a stationary point. The situation is described by a 3D version of drawing (4.2.17). Note that the intersection of the two constraint su rfaces might result in multiple curves, and each must then be considered. For example, the intersection of two thin (prolate) ellipsoids might be two closed curves. If the intersection of the constraint surfaces is null ( 1st ellipsoid inside 2nd ), the problem has no solutions. 4.4 The General Case Here we su mmarize what has been demonstrated in the above examples. We have the general gradient equation (4.1.2) emerging from the Lagrange Multiplier analysis, ∇f(r) + λ 1∇a(r) + λ2 ∇b(r) + ...... + λS-1 ∇q(r) = 0 ∇ = (∂1,∂2...∂N) . (4.1.2) This equation only exists if r is an extremum!! Solving for ∇f(r) and dotting with d r gives df( r) = ∇ f(r) • dr = – λ 1[∇a(r) • dr] – λ2 [∇b(r) • dr] – ...... – λS-1 [∇q(r) • dr] . (4.4.1) If dr at point r is chosen so all the constraints are respected (that is, we stay on all the constraint surfaces as we displace d r from r) , then ∇ a(r) • dr = 0, ∇b(r) • dr = 0 ..... and ∇ q(r) • dr = 0. Then df( r) = ∇ f(r) • dr = – λ1[0] – λ2 [0] – ...... – λS-1 [0] = 0 . (4.4.2) Since df( r) = 0, the point r is a stationary point of f(r ) subject to these constraints. The point r is then a candidate for being an extremum point of f( r). 28 5. Example 3: Extremal distances between an ellipse and a circle Many examples of the use of the Method of Lagrange Multipliers can be found in texts and on the web. See for example Trench or the wiki Lagrange Multiplier page. Our Example 3 has two constraints so S = 3. There are four coordinates so N = 4. It happens that these four coordinates are not those of a single point in E4 but represent two points in E2. This example is a bit more complex than the toy examples in texts and we obtain a numerical solution after carrying out the analytic steps of the Lagrange Multiplier algorithm. Consider then these two curves in E 2: x 2/A2 + y2/B2 = 1 ellipse centered at the origin, semimajor/minor axes A and B (x'-α)2 + (y'-β )2 = R2 . circle of radius R centered at ( α,β) (5.1) What is the minimum and maximum distance between these two curves? The distance squared is given by f = (x-x') 2 + (y-y')2 = f(x,y,x',y') N = 4 coordinates (5.2) while the constraint functions are a(x,y,x',y') = x 2/A2 + y2/B2 - 1 a(x,y,x',y') = 0 b(x,y,x',y') = (x'-α )2 + (y'-β)2 - R2 b(x,y,x',y') = 0 . (5.3) Notice that the vector r for this problem is r = (x,y,x',y'). The various derivatives are, f 1 = ∂xf = 2(x-x') a 1 = ∂xa = 2x/A2 b1 = ∂xb = 0 f2 = ∂yf = 2(y-y') a 2 = ∂ya = 2y/B2 b2 = ∂yb = 0 f3 = ∂x'f = -2(x-x') a 3 = ∂x'a = 0 b 3 = ∂x'b = 2(x'-α) f4 = ∂y'f = -2(y-y') a 4 = ∂y'a = 0 b 4 = ∂y'b = 2(y'-β ) . (5.4) The H function of (2.1) is then, H(x,y,x',y', λ 1,λ2) = f + λ1a + λ2 b = (x-x')2 + (y-y')2 + λ1 [x2/A2 + y2/B2 - 1] + λ2 [(x'-α)2 + (y'-β )2 - R2] . (5.5) Compute the four first derivatives as in (2.2), 29 H1 = f1 + λ1a1 + λ2 b1 = 2(x-x') + λ 1 2x/A2 H2 = f2 + λ1a2 + λ2 b2 = 2(y-y') + λ 1 2y/B2 H3 = f3 + λ1a3 + λ2 b3 = -2(x-x') + λ22(x'-α) H4 = f4 + λ1a4 + λ2 b4 = -2(y-y') + λ 22(y'-β) . (5.6) Set these derivatives to 0 as in (2.4) to find this set of six equations (2.6) (x-x') + λ 1 x/A2 = 0 1 (y-y') + λ 1 y/B2 = 0 2 - (x-x') + λ2 (x'-α) = 0 3 - (y-y') + λ 2 (y'-β) = 0 4 x 2/A2 + y2/B2 - 1 = 0 5 (x'-α)2 + (y'-β)2 - R2 = 0 6 ( 5 . 7 ) where the 6 unknowns are x,y,x',y',λ 1,λ2 . Notice that a possible λ i solution is λ1= λ2 = 0 in which case one gets x = x' and y = y'. Certainly this is an extremum situation since distance squared is f = (x-x')2 + (y-y')2 = 0, a minimum. Insertion of x' = x and y' = y into 5 and 6 above and elimination of x yields an equation for y of the form k1y4 + k2y3 + k3 y2 + k4 y + k4 = 0. x = ± A 1-(y/B )2 (5.8) If the ellipse and circle intersect, one will find that the above equation has 1,2,3 or 4 real roots which then represent the intersection point(s) of the two curves. If the curves don't intersect, the solution y values will be complex so there are no physical solutions for λ1= λ2 = 0 . So how does one go about solving the set of 6 equations shown in (5.7)? The possible first step is to find a viable set of λ i. From a version of (2.8) using the first and third rows of (2.7) one may write, - ⎝⎛ ⎠⎞f1 f3 = ⎝⎛ ⎠⎞ a1 b1 a3 b3 ⎝⎛ ⎠⎞λ1 λ2 ≡ M ⎝⎛ ⎠⎞λ1 λ2 (2.8) where M = ⎝⎛ ⎠⎞ a1 b1 a3 b3 = ⎝⎛ ⎠⎞ a1 0 0 b3 ( 5 . 9 ) so that M -1 = ⎝⎛ ⎠⎞ b3 0 0 a1 / det(M) = 1 a1b3 ⎝⎛ ⎠⎞ b3 0 0 a1 . ( 5 . 1 0 ) Therefore, 30 - ⎝⎛ ⎠⎞λ1 λ2 = M-1 ⎝⎛ ⎠⎞f1 f3 = 1 a1b3 ⎝⎛ ⎠⎞ b3 0 0 a1 ⎝⎛ ⎠⎞f1 f3 ( 5 . 1 1 ) and thus -λ1 = 1 a1b3 b3f1 - λ2 = 1 a1b3 a1f3 ⇒ λ1 = -f1/a1 = -2(x-x')/[2x/A2] = - A2(x-x')/x λ 2 = - f3/b3 = 2(x-x')/[2(x'-α)] = (x-x')/(x'-α ) and we then have a candidate solution set { λi}. The first four equations of (5.7) become (x-x') + [- A 2(x-x')/x] x/A2 = 0 1 (y-y') + [- A2(x-x')/x] y/B2 = 0 2 - (x-x') + [ (x-x')/(x'- α)] (x'-α) = 0 3 - (y-y') +2[ (x-x')/(x'- α)] (y'-β) = 0 4 or (x-x') - (x-x') = 0 1 (y-y') - (A/B) 2(x-x')y/x = 0 2 - (x-x') + (x-x') = 0 3 - (y-y') + (y'- β)(x-x')/(x'- α) = 0 4 . (5.12) Two of these equations are identities, while the other two are (y-y') - (A/B) 2(x-x')y/x = 0 2 - (y-y') + (y'- β)(x-x')/(x'- α) = 0 4 or (y-y')x - (A/B)2(x-x')y = 0 2 - (y-y')(x'-α ) + (y'-β)(x-x') = 0 4 . (5.13) Including the last two equations (5.7), one has four equations in four unknowns x,y,x',y' : (y-y')x - (A/B) 2(x-x')y = 0 2 - (y-y')(x'-α ) + (y'-β)(x-x') = 0 4 x 2/A2 + y2/B2 - 1 = 0 5 (x'-α)2 + (y'-β)2 - R2 = 0 . 6 (5.14) This is a system of four 2nd degree polynomial equa tions in four variables. When we ask Maple to analytically solve these equations for the four unknowns x,y,x',y' the results are very messy and involve roots of fourth degree polynomials. In this situation, the analytic solution exists, but is so complicated it 31 hardly seems very useful. We selected this example just to show that even a relatively simple problem can be nearly intractable analytically. We turn then to a numerical solution for the particular ellipse + circle case shown in this scaled figure where the background boxes are 1 unit squares, (5.15) The circle has radius R = 2 and is centered at ( α,β) = (3,1). The semi ellipse axes are A = 8 and B = 4. We enter into Maple the 6 original equations (5.7) in the 6 unknowns x,y,x',y', λ1,λ2 (xp = x', L1 = λ1, etc) : ( 5 . 1 6 ) We next set in the specific parameters for the figure shown above ( 5 . 1 7 ) 32 Without given a starting point, Maple finds the following solution which we approximate as λ 1 = -2.8 λ2 = .32 (x,y) = (3.66, 3.56) (x', y') = (3.50, 2.94) . (5.18a) Plotting this on our figure, (5.18b) we see that this solution corresponds to the minimu m distance between the two curves. Maple computes this minimum distance to be .6406470846 units. Another solution is found by giving Maple a search starting point (-8,0) on the ellipse and (5,1) on the circle : which we approximate as λ 1 = -104 λ2 = -6.5 (x,y) = (-7.99, -.22) (x', y') = (4.99, 1.22) . (5.19a) Plotting this solution, 33 (5.19b) we see that this solution corresponds to the maximu m distance between the two curves. Maple computes this maximum distance to be 13.05541338 units. A third candidate solution is found by giving Maple a starting point (-8,0) on the ellipse and (1,1) on the circle: which we approximate as λ 1 = -72 λ2 = 4.5 (x,y) = (-7.99, -.22) (x', y') = (1.01, 0.78) . (5.20a) Plotting this solution, (5.20b) one sees that this solution represents a "stationary point " but is not in fact a solution to the problem. Finally, we find this fourth solution, which we approximate as 34 λ1 = -20.2 λ2 = -2.3 (x,y) = (3.66, 3.56) (x',y') = (2.50, -.94) . (5.21a) Plotting this solution, (5.21b) we find a situation similar to the third candidate so lution -- a stationary point which is neither a global maximum nor a global minimum. 35 6. Example 4: The Boltzmann factor in Statistical Mechanics 6.1 Statement of The Boltzmann Extremum Problem I magine a system which has m state levels into which particles can be placed. At a given energy level εi, the number of available states is gi, known as the degeneracy at energy εi. A particular state of the entire system of particles is characterized by the number of particles N i in each of the states, so system state = (N 1,N2...Nm) = N , a vector. The total number of particles in the system is fixed at M = ΣiNi, and the total energy is U = ΣiεiNi. The system is isolated so its initi al energy U does not change, nor can particles be created or destroyed, so U and M are both fixed constants. The system of particles (perhaps a box of gas atoms) is assumed to be at thermal equ ilibrium so all its macroscopic characteristics (like pressure and temperature) are stable. Imagine a billion of these particle systems (an "ensem ble"). Each system settles into its own partition set {N i} = N . If M is very large, one will find that those billion N vectors are all very similar. This is so because certain {N i} partitions are statistically favored over ot her partitions just because there are a lot more "states" available to a system with those favored {N i} partitions. So, how many states are available ("accessible") to such a system characterized by vector N which represents some partitioning of the particles {N i}? The number of different "ways" of placing N indistinguishable balls into g boxes is a fascinating elementary problem and the answer is ⎝⎛ ⎠⎞ N+g-1 N , a binomial coefficient. Footnote : One arrives at this answer by thinking of N balls and g-1 "partitions" between boxes, all of which have to be laid out in a row, such as | * * * | * | | * for the case N = 5 and g-1 = 4. Here the leftmost bin is empty, as is the second bin from the ri ght. Each layout of the 9 objects is characterized by picking a committee of 5 balls from the 9 positions (or a committee of 4 partitions). So in this case the number of unique layouts is ⎝⎛ ⎠⎞ 9 5 = 9! 5!4! = ⎝⎛ ⎠⎞ 9 4 . (6.1.1) The count ⎝⎛ ⎠⎞ N+g-1 N allows multiple balls to be in any given box. This count is then associated with so- called Bose-Einstein statistics -- multiple bosons are allowed in the same state (for example, photons or any integral-spin particles). If at most one ball is allo wed in any given box , the count of ways is just the number of ways to pick a committee of N boxes from the set of g boxes, and this count is ⎝⎛ ⎠⎞ g N . This count is then associated with so-called Fermi-Dirac statistics -- only one fermion is allowed in a state (for example electrons or any half-integral spin partic les). The fact that electrons are fermions is the reason the periodic table of elements exists. However, we shall assume that g i >> Ni (known as the Boltzmann limit), and this allows a simplification of the state count. Both the Bose and Fermi state counts reduce to the same count in this limit: 36 ways = ⎝⎛ ⎠⎞ N+g - 1 N ≈ ⎝⎛ ⎠⎞ g N = g! N! (g-N)! = g(g-1)..(g-N+1) N! ≈ gN N! . (6.1.2) So in our particle problem, the number of ways of having N 1 particles in states with energy ε1 is g1N N1! . For each of these ways, there are g2N N2! ways of putting N 2 particles in energy level ε2. And so on. Thus, the total number of available states ("microstates") for a partition {N i} of the M particles is given by Ω(N1,N2...Nm) = g1N1 N1! g2N2 N2! ....gmNm Nm! . (6.1.3) So here is our initial extremum problem : Find N which maximizes Ω(N) = g1N1 N1! g2N2 N2! ....gmNm Nm! subject to these two constraints: Σ iNi = M and ΣiεiNi = U . (6.1.4) When M is a very large number (like Avogadro's number), it turns out that Ω(N) has a very strong and sharp maximum at a certain N which is the solution value N for this extremum problem. Thus, when one examines an ensemble of such systems, this solution N is the overwhelmingly likely partitioning of the particles into {N i}. Our task is then to solve this problem for N. This vector N = (N1, N2...Nm) is the vector r = (x1,x2...xN) of our general Lagrange Multiplier presentation in Section 2, so in this application the nu mber of components in the vector is N = m, and we don't want to confuse this previous use of symbol N with N or Ni of the present problem! One might observe that the x i in the general analysis were reals, whereas the N i here are integers. We can dispense with this issue by simply writing N i! = Γ (Ni+1) which then serves to interpolate N i between its integer values. Another approach is to replace N i everywhere in the above problem statement by ni ≡ Ni/M. Since M and the solution N i values are very large integers, one can regard this n i as essentially a continuous real (albeit rational) variable. It should be noted that the value of U is restricted to a certain range. Mε min ≤ U ≤ Mεmax . ( 6 . 1 . 5 ) The limits represent the two extreme possible partitions {N i} of the particles (all in lowest energy state, or all in highest energy state). It is useful to define the average energy of a particle, u ≡ U/M = average energy of particle in the system (6.1.6) and then (6.1.5) says, ε min ≤ u ≤ εmax . ( 6 . 1 . 7 ) 37 6.2 Solution of The Boltzmann Extremum Problem Maxim izing the state count Ω of (6.2) is the same as maximizing f ≡ lnΩ, and from (6.1.3), f(N 1, N2...Nm) ≡ lnΩ(N1,N2...Nm) = Σi [ Ni ln gi ] - Σi [ ln Ni! ] (6.2.1) where Σi means Σi=1m. Assuming all the N i are very large numbers, we may use Stirling's formula, x! ≈ 2πx xx e-x ⇒ lnx! ≈ ln 2π + (1/2) lnx + xlnx - x ≈ xlnx - x, x >> 1 . (6.2.2) Then ln N i! ≈ NilnNi- Ni so that, f = ln Ω = Σi [ Ni lngi ] - Σi [ ln Ni! ] ≈ Σi [ Ni lngi ] - Σi [NilnNi- Ni] = Σi [(1+lngi)Ni - NilnNi ] = Σi Ni [1 + ln(g i/Ni)] . (6.2.3) Here then is a restatement of the extremum problem given in (6.1.4): Find N which maximizes f( N) = Σ i [(1+ln g i)Ni - NilnNi ] subject to these two constraints: a ( N) = ΣiNi - M = 0 b ( N) = ΣiεiNi-U = 0 . (6.2.4) Following now the prescription of Section 2, we write our "Lagrangian" H as H( N,λ) ≡ f(N) + λ 1a(N) + λ2 b(N) = Σ i [(1+ln g i)Ni - NilnNi ] + λ1 [ ΣiNi - M] + λ2 [ΣiεiNi-U] . (6.2.5) The next step (b) is to compute the m derivatives with respect to N i : f i = ∂f/∂Ni = (1+lng i) - Ni (1/Ni) - lnNi = lngi - lnNi = ln(gi/Ni) a i = ∂a/∂Ni = ∂i [ ΣjNj - M] = 1 b i = ∂b/∂Ni = ∂i [ΣjεjNj-U] = εi . (6.2.6) In step (c) we set the partials of H to zero, so 0 = H i(N,λ) = fi + λ1ai + λ2 bi = ln(g i/Ni) + λ1 + λ2εi i = 1,2...m or ln(g i/Ni) + λ1 + λ2εi = 0 i = 1,2...m or 38 ln(N i/gi) = λ1 + λ2εi . i = 1,2...m (6.2.7) Step (d) just displays the set of equations one needs to solve : ln(N i/gi) = λ1 + λ2εi i = 1,2...m Σ iNi = M / / a ( N) = 0 Σ iεiNi = U / / b ( N) = 0 . (6.2.8) This is a system of m+2 equati ons in m+2 unknowns which are the N i, λ1 and λ2. To solve these equations, we first solve the first set of m equations for N i, leaving λ1 and λ2 as unknowns: ln(N i/gi) = λ1+ λ2εi or Ni = gi eλ1 eλ2εi . ( 6 . 2 . 9 ) Already a major result has emerged from our Lagrange Multiplier analysis! The population N i of states at energy level εi depends exponentially on εi. We intuitively expect higher energy levels to be less populated than lower ones, so we expect λ2 to be negative. For this reas on, and following tradition, we set λ2 = - β where β > 0. And while we're at it, we can set eλ1 = A to simplify notation. So β ≡ - λ 2 A ≡ eλ1 ( 6 . 2 . 1 0 ) and then (6.2.9) reads, N i = A gi e-βεi . ( 6 . 2 . 1 1 ) A study of the statistical mechanics of th ermal equilibrium in fact shows that β = 1 / k T ( 6 . 2 . 1 2 ) where T is absolute temperature and k is Boltzmann's constant. The factor e -βεi in (6.2.11) is referred to as the Boltzmann factor . The quantity S = k ln Ω = kf is the entropy of a system. The idea that β = 1/kT is really the definition of absolute temperature, but all these matters are not of immediate interest, so we continue onward. From (6.2.11) the probability of a particle having energy εi is given by pi = Ni M = A M gi e-βεi . ( 6 . 2 . 1 3 ) 39 Since Σipi = 1 one finds that 1 = A M [Σigi e-βεi] which then determines the derived Lagrange multiplier constant A (expressed in terms of β), A = M Σigi e-βεi . ( 6 . 2 . 1 4 ) The denominator in (6.2.14) is called the partition function for the system. Using (6.2.11) for N i, the two constraint equations shown in (6.2.4) [ ΣiNi= M and ΣiNiεi= U] read A Σigie-βεi = M A Σiεig e-βεi = U . (6.2.15) Dividing these two equations and using (6.1.6) that u = U/M, one gets this self-consistent result, u = U M = Σigiεie-βεi Σigie-βεi = (A/M)Σigiεie-βεi (A/M)Σigie-βεi = Σipiεi 1 = Σipiεi = <εi> . (6.2.16) To find the derived Lagrange multiplier constant β = -λ2, write out the second equation of (6.2.15), Σigiεi e-βεi = U/A = (U/M)(M/A) = (u) ( Σigi e-βεi) (6.2.17) so this is the equation one must solve for β. Recall that εi, gi, and u are all known quantities. For the lowest energy states one will have εi < u ( recall that u = < εj> ) so [εi - u] will be negative. For high energy states [ εi - u ] will be positive, so that is why (6.2.17) at least stands a chance of having a solution for β. As noted in (6.2.12), the solution for β indicates the temperature of the system of particles. Changing u changes the solution β , so one concludes that the mean particle energy u is a function of temperature. For general values of ε i (6.2.17) is a transcendental equation which must be solved numerically for β. If it happens that all the εi are ratios of integers, the equation can be written as a polynomial equation, but since those integers are likely to be large, the orde r of this polynomial is large and again one must do a numerical solution. Notice that the solution to our particle system problem has the same general form regardless of the number m of energy levels ε i. If we take m very large and make the energy levels εi be very closely spaced, they approach a continuum of energy values and the conclusions apply to a classical system. For small finite m, the discrete energy levels indicate a quantum mechanical system. To summarize, we have used the method of Lagrange multipliers to solve the Boltzmann extremum problem stated in either (6.1.4) or (6.2.4). Since th e problem has two constraints, there are two Lagrange multipliers λ1 and λ2 which we replaced with the numbers A ≡ eλ1 and β ≡ - λ2 . We showed how β is determined by numerically solving (6.2.17), and then A is determined by (6.2.14). Then the solution (extremum) vector N = {Ni} is given by (6.2.11) which says N i = A gi e-βεi . 40 What we have not shown is that this extremum solu tion is a maximum. nor have we shown the fact that Ω(N) has a very sharp peak at the solution value N , which then statistically forces a system to assume a value N very close to this solution N. We shall deal with these issues in the following section. 6.3 More details of the solution Start with f = lnΩ = Σi [(1+lngi)Ni - NilnNi ] = Σi Ni [1 + ln(g i/Ni)] . (6.2.3) The constraints Σ iNi= M and ΣiNiεi= U can be regarded as two equation in the two unknowns N 1 and N2 which are easily solved to obtain, N1 = + (ε2- ε1)-1[ ε2M - U + Σ i=3m (εi- ε2)Ni] = N 1(N3,... Nm) N2 = - (ε2- ε1)-1[ ε1M - U + Σ i=3m (εi - ε1)Ni] = N 2(N3,... Nm) . (6.3.1) Notice that, for i = 3,4...m, ∂N 1 ∂Ni = + (ε2- ε1)-1 [ (εi- ε2) ] = εi-ε2 ε2-ε1 ∂N2 ∂Ni = - (ε2- ε1)-1 [ (εi- ε1) ] = - εi-ε1 ε2-ε1 . (6.3.2) Now think of f as function of N 3,,,,Nm (analogous to the first line of (3.4)) , f(N 1(N3,... Nm), N2(N3,... Nm), N3, .....Nm) . (6.3.3) Compute the total derivative df/dN i for i = 3,4...n: df dNi = f1 ∂N1 ∂Ni + f2 ∂N2 ∂Ni + fi where f i ≡ ∂f ∂Ni = f 1 εi-ε2 ε2-ε1 - f2 εi-ε1 ε2-ε1 + fi (6.3.2) = ln(g1/N1) εi-ε2 ε2-ε1 - ln(g2/N2) εi-ε1 ε2-ε1 + ln(g i/Ni) . (6.2.6) (6.3.4) Installing the solution values ln(N i/gi) = λ1 + λ2εi shown in (6.2.8) we continue, df dNi = (λ1 + λ2ε1) εi-ε2 ε2-ε1 - (λ1 + λ2ε2) εi-ε1 ε2-ε1 + (λ1 + λ2εi) = λ1 [ εi-ε2 ε2-ε1 - εi-ε1 ε2-ε1 + 1] + λ2 [ ε1 εi-ε2 ε2-ε1 - ε2 εi-ε1 ε2-ε1 + εi ] 41 = λ1 [ εi-ε2 ε2-ε1 - εi-ε1 ε2-ε1 + ε2-ε1 ε2-ε1 ] + λ2 [ ε1 εi-ε2 ε2-ε1 - ε2 εi-ε1 ε2-ε1 + εi ε2-ε1 ε2-ε1 ] = λ1 [0] + λ2[0] = 0 . (6.3.5) This is the expected result since th e solution is supposed to be a critical point of f given in (6.3.3). Our main interest is in the second derivative d2f/dNi2 . If this comes out negative at the solution point, we know we have a maximum. Applying d/dN i to the second line in (6.3.4), and then using the first line of (6.3.4) sequentially for f → f1, f2 and fi , one gets d2f d2Ni = df1 dNi εi-ε2 ε2-ε1 - df2 dNi εi-ε1 ε2-ε1 + dfi dNi = [ f 11 εi-ε2 ε2-ε1 - f12 εi-ε1 ε2-ε1 + f1i] εi-ε2 ε2-ε1 f ij ≡ ∂2f ∂Ni∂Nj = fji - [ f 21 εi-ε2 ε2-ε1 - f22 εi-ε1 ε2-ε1 + f2i] εi-ε1 ε2-ε1 + [ f i1 εi-ε2 ε2-ε1 - fi2 εi-ε1 ε2-ε1 + fii] . i = 3,4....m (6.3.6) Recall now that f = ln Ω = Σi [(1+lngi)Ni - NilnNi ] (6.2.3) f i = lngi - lnNi (6.2.6) fij = - δij(1/Ni) . ( 6 . 3 . 7 ) Therefore (6.3.6) simplifies to d 2f d2Ni = f11 [ εi-ε2 ε2-ε1 ]2 + f22 [ εi-ε1 ε2-ε1 ]2 + fii = – { 1 N1 [ εi-ε2 ε2-ε1 ]2 + 1 N2 [ εi-ε1 ε2-ε1 ]2 + 1 Ni } . (6.3.8) This says that the curvature of f in all independent directions N i for i = 3,4...m is negative and this is true for all vectors N = {Ni}. This is of course then true for the solution vector N. But for the solution vector N we also have df/dN i = 0 from (6.3.5) and ther efore the solution is a maximum of f and therefore of Ω. Finally, we may estimate the width of the peak of f and then of Ω. The energy fractions in (6.3.8) are on the ballpark order of unity and the N i are on the order of M, the total number of particles. If we do a Taylor expansion of f( N) in the i direction only (i = 3,4...m) about the solution N, then 42 f(Ni + dNi) ≈ f(Ni) + (df/dN i)dNi + (1/2)(d2f/d2Ni) (dNi)2 . (6.3.9) But at the solution point (df/dN i) = 0 so one finds df = f(N i + dNi) - f(Ni) ≈ (1/2)(d2f/d2Ni) (dNi)2 ( 6 . 3 . 1 0 ) so (dN i)2 ≈ 2 |df| |d2f/d2Ni| . ( 6 . 3 . 1 1 ) Suppose we are interested in finding dN i such that the Ω(N) drops to half its peak value. Then Ωhalf/Ωpeak = 1/2 ⇒ ln Ωhalf - ln Ωpeak = 1/2 ⇒ fhalf - fpeak = ln(1/2) ⇒ |df| = f peak - fhalf = ln(2) = 0.7 . (6.3.12) In this case, the half-width (6.3.11) of the Ω plot is dNi = 1.4 |d2f/d2Ni| . ( 6 . 3 . 1 3 ) For a ballpark estimate, we can interpret (6.3.8) as saying d2f d2Ni ~ -3/M ⇒ dN i ≈ 1.4 M 3 ~ M (6.3.14) so the fractional half-width is roughly. dN i M ~ 1 M . ( 6 . 3 . 1 5 ) For a macroscopic box of gas particles, one might have M = 1020 and then dNi M ~ 10-10 . Since this is the fractional half-width of the peak in Ω, the peak is extremely sharp. This then is why, in an ensemble of macroscopic particle systems, all systems will have almost the same N value. 43 6.4 A numerical example with plots For this example we choose three energ y levels (m = 3) so the state count is given by Ω(N1,N2,N3) = g1N1 N1! g2N2 N2! g3N3 N3! (6.1.3) (6.4.1) and correspondingly (this uses the Stirling approximation), f = ln Ω = Σi=13 Ni [1 + ln(g i/Ni)] . (6.2.3) (6.4.2) Assume the three energy levels are ordered in this manner, where ε 1 is the lowest energy state, ε 1 < ε2 < ε3 . and recall u = U/M and ε1 ≤ u ≤ ε3 . (6.4.3) We regard both N 1 and N2 as functions of N 3 where, according to (6.3.1), N1 = + (ε2- ε1)-1[ ε2M - U + ( ε3- ε2)N3] = N 1(N3) N2 = - (ε2- ε1)-1[ ε1M - U + ( ε3 - ε1)N3] = N 2(N3) . (6.3.1) (6.4.4) It seems clear that one must have 0 ≤ N 1 ≤ M 0 ≤ N2 ≤ M . ( 6 . 4 . 5 ) Using the expressions (6.4.4) in these inequalities and making use of (6.4.3) and the fact that ε 1 ≤u ≤ ε3, one finds that all four inequalities can be summarized as just two inequalities which we write as, N3min ≤ N3 ≤ N3max N 3min = max(u-ε2 ε3- ε2 M, 0) N 3max = min( u-ε1 ε3- ε1 M, M) = u-ε1 ε3- ε1 M . // since u ≤ ε3 (6.4.6) To find the Lagrange multiplier β we must solve (6.2.17), Σi=13 [ εi - u ] gi e-βεi = 0 . (6.2.17) (6.4.7) Then the other derived multiplier is given by (6.2.14) A = M Σ igi e-βεi = M Z where Z ≡ Σigie-βεi . (6.2.14) (6.4.8) The components N i of the solution vector are then given by 44 Ni = A gi e-βεi . (6.2.11) (6.4.9) The curvature of f at its peak is, d 2f d2N3 = – { 1 N1 [ ε3-ε2 ε2-ε1 ]2 + 1 N2 [ ε3-ε1 ε2-ε1 ]2 + 1 N3 } . (6.3.8) (6.4.10) The half-width of the Ω peak is then given by dNi = 1.4 |d2f/d2Ni| . (6.3.13) (6.4.11) For our specific Maple plots, we shall assume ε3= 3 M = 1000 = number of particles ε2 =2 u = 1.5 = average energy ε1= 1 g 1= g2= g3 = 5000 = degeneracy (6.4.12) Here then is our Maple code. We first enter Ω, The quantities N 1 and N2 are then replaced as in (6.4.4) [ we use e i in place of εi] , We next enter most of the data shown above and use it to compute N 3min and N3max 45 so the legal range of N 3 is (0,250). Here is a plot of N 1 (red) and N 2(black) versus N 3 : (6.4.13) We then enter the degeneracies and take a look at Ω in its numerical form, (6.4.14) As is typical in statistical mechanics, at its peak the value of Ω is a very large number. Even with the relatively small number of particles M = 1000 and dege neracy g = 5000, the peak value is (see below) Ωpeak = 6.71 x 101519 ( 6 . 4 . 1 5 ) Maple is uncomfortable plotting numbers larger than about 10 40 so we pre-scale Ω down by 101519 to make the scaled Ω be Maple-digestible. This scaling process in turn creates numbers smaller than 10-40 which are also rejected by the Maple plotter, so we use a Heaviside function to pin small numbers to 0 unless they are greater than the arbitrary value .01. Here then is a plot of Ω(N3), 46 ( 6 . 4 . 1 6 ) We broke the display command into two pieces: first create xx then plot xx. Replacing : with ; after the plot command gives one a list of numbers to be plo tted, and from this list one can determine a reasonable scale factor. Once that is found, the plotting c oordinates need not be displayed. The curve Ω(N 3) has the typical Gaussian (normal) shape which no doubt the en ergetic reader could show results from the central limit theorem. A plot of f = ln( Ω) has a much smoother and broader shape, ( 6 . 4 . 1 7 ) Recall from (6.3.8) that d2f d2Ni is everywhere negative, consistent w ith the cupping down seen in this plot. Next, we have Maple solve (6.2.17) for the Lagrange multiplier β, 47 With w = e-β the above equation is 3w2+w-1 = 0 which has solutions w = ( -1 ± 13)/6. The minus sign gives the complex β noted, while the plus sign gives the physical solution β = - ln[( -1 + 13)/6] = .834. From this β value we compute the other Lagrange parameter A from (6.2.14), Once these parameters are known, we use (6 .2.11) to display the solution vector N = {Ni}. Since we set u < ε 2, we see here a normal distribution where high er energy states have fewer particles. Our final task is to use (6.3.8) and (6.3 .13) to determine the half-width of the Ω peak, which says Δ N1/2 ≈ 7.5. This width and the solution N 3 = 116.2 can be verified from this blowup plot of the Ω peak, 48 ( 6 . 4 . 1 8 ) Notice that this Maple code does not use the Stirling approximation and provides a continuous interpolation of N! as Γ(N+1). The peak value of Ω is found to be In the above example we used u = 1.5 < ε2 = 2. If we instead use u = 2.5 > ε2 = 2, we find that the state populations are as shown above but in reverse order. Such a inverted situation corresponds to a negative absolute temperature T ( β = -.834 = 1/(kT) ) and cannot therefor e arise in thermal equilibrium. It does arise in a laser medium which is driven by some power source ("population inversion"). 49 Appendix A : Matrix rank equals the number of independent columns and rows This Appendi x proves that, for a general m x n matrix M, the number of linearly independent columns and the number of linearly independent rows are the same, and both numbers are equal to the rank of the matrix. The reader will appreciate that this fact is non-obvious. Shilov proves it but the proof is a bit spread out over several sections. Here we provide a self-contained alternate proof. We first quote a series of square-matrix theorems th at are proved from scratch in our Appendix B. These are doubtless very familiar to the reader, but we just want to get them stated. The derivations of these theorems may be less familiar. Theorem 1 : det(M T) = det(M). Switching rows with column s does not change a determinant. (B.1.10) Theorem 2: det(M) can be represented in these two ways, where ε is the permutation tensor (B.2.1) : det(M) = Σ a1a2... an εa1a2... an M1a1M2a2 ...... Mnan (B.2.10a) det(M) = Σ a1a2... an εa1a2... an Ma11Ma22 ...... Mann . (B.2.10b) Theorem 3: For a square matrix M, adding a multiple of one row to another does not change det(M). The same is true for adding a multiple of one column to another. (B.2.12) Theorem 4: Swapping two rows (columns) of square matrix M causes det(M) → - det(M). (B.2.13) Corollary 4: If two rows (columns) of square matrix M are the same, then det(M) = 0. (B.2.14) Theorem 5 : det(AB) = det(A)det(B) (B.2.15) Corollary 5 : det(ABC) = det(A)det(B)det(C) and so on. (B.2.16) Theorem 6 (Cofactor Expansions): det(M) = Σ n Msn cof(Msn) = Σn (rs)n cof(Msn) work across row s s = 1,2....n (B.3.16a) det(M) = Σn Mns cof(Mns) = Σn (cs)n cof(Mns) work down column s s = 1,2....n (B.3.16b) Theorem 7: M-1 = CT det(M) = [cof(M)]T det(M) = cof(MT) det(M) for square matrix M. (B.4.8) Theorem 8 (Cramer's Rule): y = M x ⇒ xs = det(M[ cs→ y]) det(M) (B.4.14) To this list, we shall now add four new Theorems which will be proven right here. The definition of linear independence and linear dependence is given in the discussion surrounding (1.6) and won't be repeated here. 50 Theorem 9: Columns of square matrix M are linearly dependent ⇔ det(M) = 0 (A.1) Contrapositives : Columns of square matrix M are linearly independent ⇔ det(M) ≠ 0 If we can prove Theorem 9 as stated, then the theorem is also true fo r rows. The reason is that swapping rows and columns corresponds to M ↔ MT and Theorem 1 says det(MT) = det(M). Proof of ⇒: According to (1.6) and nearby discussion, our premise of linear dependence says that either one or more columns of M are zero, or at least one column can be written as a linear combination of the other columns, so perhaps cj = Σikici . If a column is zero, det(M) = 0 from Theorem 6 going down that column. If cj = Σi≠jkici, one can add - Σi≠jkici to cj without changing det(M) according to Theorem 3 . But this makes the new column j vanish, so again det(M) = 0. Proof of ⇐ : Here we shall prove the contrapositive instead of the claim: claim: det(M) = 0 ⇒ the columns of M are linearly dependent contrapositive: the columns of M are linearly independent ⇒ det(M) ≠ 0 If the column vectors c i of matrix M are linearly independent, then (1.6) says Σ xici = 0 ⇒ x = 0. Therefore we know that ΣxiMji = 0 ⇒ xj = 0 which is the same as Mx = 0 ⇒ x = 0. Thinking of M:En → En, we claim that mapping M x = y is one-to-one. Certainly each x goes into a single y, Could a y map back into two different values of x, call them x and x' with x ≠ x' ? Then M x = y and M x' = y . Subtract to get M( x-x') = 0. But since M z = 0 ⇒ z = 0 (as just shown), we find x = x' . Thus, the mapping Mx = y really is one-to-one, and that means it is invertible and the inverse is unique. From Theorem 7 the inverse is in fact given by M-1 = cof(MT)/det(M). For M-1 to exist, one must have det(M) ≠ 0. Another proof is very simple but perhaps less c onvincing. One can show (see our Tensor Analysis document (B.6.10) in Refs) that if a set of column vectors ci of matrix M spans an n-piped in En, the "volume" of that n-piped is given by V = |det(M)| . If those vectors are linearly independent, then V ≠ 0. This is obvious in 2D (volume = area) and 3D but less obvious for general En. For example, in 3D if a third vector lies in the plane of the other two (is dependent), there is no volume. Theorem 10 (mini-columns theorem). Let R be a subset of N rows of nxm matrix A, and let S be a subset of columns {c i} in A. The matrix elements of A included in the intersection of these two sets form a set of "mini-columns" { Ci}. This set of mini-columns forms a matrix B within A. Matrix B could be contiguous, or it could be non-contig uous in one or both directions. Th e claim and its contrapositive are: (a) { Ci} linearly independent in B ⇒ {ci} linearly independent in A (b) {c i} linearly dependent in A ⇒ {Ci} linearly dependent in B (A.2) A picture is worth a thousand words. Here matrix B happens to be contiguous in both directions. 51 (A.3) Proof of (b): If the set { ci} is linearly dependent in A, then from (1.6) one can write Σi∈S kici = 0 with k ≠ 0. This is a set of m equations which contains as a subset the set of N equations Σi∈S kiCi = 0. This last equation then says that the { Ci} are linearly dependent in B. There is of course a similar "mini rows theorem". Theorem 11 . If all kxk minors in k columns of matrix A vanish, the k columns are linearly dependent. ( A . 4 ) Contrapositive : If k columns are linearly independent, they must contain at least one non-zero kxk minor. Once proven for columns, replacing A →AT gives the same theorem for rows. Proof : Assume matrix A has m rows and n columns. In order for kxk minors to exist, one must have k ≤ min(m,n). Gather up the k columns and make them be the leftmost k columns of a new matrix B. Since these columns have m elements, and since k ≤ m, add m-k arbitrary new columns to the right of the k columns so that matrix B is then a square matrix. On the left below we show the k columns taken from matrix A in red, and then the arbitrary added columns are shown in blue. (A.5) Consider now the process of computing det(B) by going down the rightmost column using the standard cofactor sum formula of Theorem 6 . This det(B) is a linear combinati on of (m-1)x(m-1) minors all in the left m-1 columns. Consider one of these minors. If we evaluate it using the same cofactor formula (of one less dimension), it will be a linear combination of (m-2)x(m-2) minors in the left m-2 columns. We keep going until we are evaluating a set of kxk cofactors in the left k columns. But these all vanish by the theorem premise. Thus, reversing this logic we conclude that det(B) = 0. The picture on the right above shows one minor at each level of the descent just described (red dot goes with red minor, etc). Now suppose it were possible that the k red columns were independent. Since the added blue columns are arbitrary, we could certainly find a set of added blue columns so that all m columns were independent. 52 But then we would have det(B) ≠ 0 by Theorem 9 . But we just showed that det(B) = 0, so it must not be possible to have the k red columns be i ndependent, so they must be dependent. We are finally ready for the Main Act. Recall first that : Definition : The rank r of an m x n matrix A is the dimension of the largest non-vanishing minor within A. [Shilov 1.92 ] Thus, r ≤ min(m,n). (1.5) (A.6) Theorem 12 . For a general nxm matrix A, (A.7) rank(A) = r ⇔ A has exactly r linearly independent columns Once this theorem is proved for columns, it is also true for rows since rank(A T) = rank(A). Proof of ⇒ : If rank(A) = r, A must have at least one non-vanishing rxr minor. Think of this minor being a set of mini-columns as in Theorem 10. Since then det(minor) ≠ 0, this set of mini-columns is linearly independent according to Theorem 9 . Thus, the set of corresponding full columns (those that pass down through the minor) are linearly independent from Theorem 10 (a). So matrix A has at least these r independent columns. Now let k = r+1. Since rank(A) = r we know that all kxk minors in A vanish. By Theorem 11, any set of k columns must be linearly dependent. Since k = r+1, any r+1 columns are linearly dependent, so the number of independent columns is just r, the number of columns passing down through the minor. Proof of ⇐ : If there are r linearly independent columns, th en at least one rxr minor within those columns must be non-zero (contrapositive of Theorem 11). Suppose there were an (r+k)x(r+k) non-vanishing minor. Then the r+k mini-columns of that minor would be independent, and thus by Theorem 10 there would be r+k independent full columns. But the th eorem premise says there are only r independent columns, so all minors larger than rxr must vanish. Therefore rank(A) = r. Sometimes the number of linearly independent rows of a matrix is called the row rank , and the number of independent columns is called the column rank . Theorem 12 and its row version show that column rank = row rank = rank = dimension of largest non-vanishing minor (A.8) Basis Minors and Basis Columns. If rank(A) = r, we know there must exist at least one non-vanishing rxr minor in A. Shilov refers to such a minor as a basis minor and the columns passing through this minor are called basis columns. We showed in Theorem 10 that the set of such basis columns is linearly independent (which is why they are called basis columns). Obviously any of these columns can be written as a linear combination of th e basis columns, such as c2 = Σi=1k kici = 1 c2. If k = r+1, all kxk minors vanish and by Theorem 11 all sets of k = r+1 columns are linearly dependent. Thus, every non-basis column is a linear combination of the basis columns. We have thus proved Shilov's Basis Minor Theorem which we quote from section 1.93 of his book (p 25) (A.9) In his proof that column rank = rank, Shilov uses the above theorem as a starting point. 53 Appendix B : Determinants There must be hundreds of definitions and theorems inhabiting Matrix World. Here we provide our own "take" on a small number of them relating to the determinant of a square matrix. The reader might be surprised to see "group theory" playing a major role in the development below. There are many ways to skin a cat, and the following thread is just one of those ways. Most of the Theorems proven below are used in Appendix A. B.1 Definition of the determinant First, we define z 0 to be the n-component vector of increasing integers 1 to n, z0 ≡ ⎝⎜⎛ ⎠⎟⎞ 1 2 ... n . ( B . 1 . 1 ) Let a be some permutation (reordering) of these integers, so write a ≡ ⎝⎜⎛ ⎠⎟⎞ a1 a2 ... an = A ⎝⎜⎛ ⎠⎟⎞ 1 2 ... n = A z 0 ( B . 1 . 2 ) where A is an nxn matrix which has 1's in th e right places to create this permutation vector a. For example, a = Az o ↔ ⎝⎜⎛ ⎠⎟⎞ 1 3 2 = ⎝⎜⎛ ⎠⎟⎞ 1 0 0 0 0 1 0 1 0 ⎝⎜⎛ ⎠⎟⎞ 1 2 3 . The vector a can be obtained from the vector z 0 by making S a pairwise swaps of elements of z0. Although S a is not unique, the number (-1)Sa is unique. For example, to get from (1,2,3) to (1,3,2) one could swap the second pair so S a = 1, but one could then swap the first pair twice and then S a = 3. This number (-1)Sa is of course ±1 and we shall call it the parity of the permutation a, Parity( a) ≡ (-1)Sa . ( B . 1 . 3 ) Fact: Parity(C a) = Parity(CA z 0) = (-1)Sc+Sa = Parity(C-1a) (B.1.4) Proof : Doing permutation CA involves first doing A with its S a swaps, and then doing C with its S c swaps, for a total of S a+Sc swaps. The inverse permutation C-1 obviously involves the same number of swaps as the permutation C. 54 Moving now toward the definition of the determinant of a matrix M, define ( Π meaning "product") Π(a,b; M) ≡ Mab ≡ Ma1b1Ma2b2 ...... Manbn = product of n factors . (B.1.5) The notation Π(a,b; M) is easier to deal with than M ab, but we shall use both these notations. Suppose c = Cz 0 is some arbitrary permutation of z0. Then: Fact: Π(Ca,Cb; M) = Π(a,b; M ) ( B . 1 . 6 ) Proof : Applying the same permutation to both the a k and bk indices of M a1b1Ma2b2 ...... Manbn, just reorders the terms in the product but the product stays the same. Fact : Π(a,b; MT) = Π(b,a; M ) ( B . 1 . 7 ) Proof : Π(a,b; MT) = MT a1b1MT a2b2 .... MT anbn = Mb1a1Mb2a2 .... Mbnan = Π(b,a; M) We shall now define the determinant of an n x n matrix M in the following admittedly obscure manner (later we will show that it reduces to more familiar forms), det(M) ≡ 1 n! Σa Σb (-1)Sa+Sb Π(a,b; M) = 1 n! Σa,b (-1)Sa+Sb Mab = 1 n! Σa Σb (-1)Sa+Sb Ma1b1Ma2b2 ...... Manbn . (B.1.8) Here Σa means the sum over all permutations a of z0. In (B.1.8) the columns and rows of M are on a completely equal footing. Notice that a and b are in effect dummy summation indices. If we do a↔b, the expression for det(M) is unchanged. Thus one can rewrite (B.1.8) as det(M) = 1 n! Σb Σa (-1)Sb+Sa Π(b,a; M) . (B.1.9) We are now ready for our first determinant theorem: Theorem 1 : det(M T) = det(M). Switching rows with columns does not change a determinant. (B.1.10) Proof : det(MT) = 1 n! Σb Σa (-1)Sb+Sa Π(b,a; MT) // (B.1.9) applied to MT = 1 n! Σb Σa (-1)Sb+Sa Π(a,b; M) // (B.1.7) applied to MT = det(M) // (B.1.8) 55 B.2 The permutation group and the permutation tensor ε Definition. A set of ele ments {g i} form a group G if : • gigj is also in the group (closure) • gi-1 exists for each g i in G, where g i-1 is also in G (inverse) • (gigj)gk = gi(gjgk) (associative) (B.2.1) Comment : It is implicit in the above definition that th e group has some "operation" which gives meaning to gigj, which we shall just think of as "multiplica tion". When group elements are represented by matrices, that operation is multiplication of those ma trices, and that will apply to the permutation group below. Fact : giG = G (the rearrangement theorem) (B.2.2) This says that multiplication of all the elements of a group by an element g i in the group creates a reordering of the group elements. Proof . Consider g iG = gi [g1, g2....gn] = [gig1, gig2....gign] = set of n elements. Unless two elements are the same, this must exhaust the entire group. How do we know that g ig1 and gig2 might not be the same? Apply g i-1 from the left and that would say g 1 = g2 which is not the case. Fact: Σ g f(g) = Σ g f(g1g) if Σg runs over the entire group G (B.2.3) Proof: In Σg f(g1g), as g runs over G, the argument g' ≡ g1g runs over G by the rearrangement theorem. Thus, the sum Σ g f(g1g) is just a reordering of the terms in the sum Σg f(g). It is easy to show that the set of permutations of z 0 forms a group and therefore the above facts can be used. For example, the product of two permutations is a permutation, and ever y permutation clearly has an inverse, and (AB)C = A(BC) for any matrices. Fact (B.2.3) can be written ΣB f(Bzo) = ΣB f(CB z0) // b = Bz0 Here the sum Σ B is over all permutation matrices which correspond to permutations b of z0, and C is some arbitrary permutation matrix associated with some permutation c of z0. An equivalent way of stating the above involves a direct summation Σb over all permutation vectors b of z0 , Σ b f(b) = Σb f(Cb) . // permutation sum rearrangement theorem (B.2.4) Throwing in the arbitrary permutation C merely causes a reordering of the sum. This is an extremely powerful and useful fact. Consider then : 56 det(M) ≡ 1 n! Σa [Σb (-1)Sa+Sb Π(a,b; M) ] // (B.1.8) = 1 n! Σa [Σb (-1)Sa+Sb Π(A-1a,A-1b; M)] // (B.1.6) with C = A-1; A-1a = A-1Az0 = z0 = 1 n! Σa [ Σb Parity(A-1b) Π(z0,A-1b; M) ] // (B.1.4) with C = A-1 and a = b = 1 n! Σa [ Σb Parity( b) Π(z0,b; M) ] // (B.2.4) with C = A-1 (key step) = Σb Parity( b) Π(z0,b; M) // n! identical terms in Σa = Σa Parity( a) Π(z0,a; M) // rename dummy sum variable = Σa Parity( a) M1a1M2a2 ...... Mnan . // (B.1.5) (B.2.5) Since det(MT) = det(M) from Theorem 1 (B.1.10), we can also write this result as det(M) = Σa Parity( a) Ma11Ma22 ...... Mann ( B . 2 . 6 ) Proof : det(M) = det(MT) = Σa Parity( a) MT 1a1MT 2a2 ...... MT nan = Σa Parity( a) Ma11Ma22 ...... Mann Here then are last two results: ( recall that S a is the number of pairwise swaps to get from z0 to a ) det(M) ≡ Σa Parity( a) M1a1M2a2 ...... Mnan Parity( a) = (-1)Sa (B.2.7a) det(M) ≡ Σa Parity( a) Ma11Ma22 ...... Mann Parity( a) = (-1)Sa (B.2.7b) Definition : The permutation tensor εijk.. (n subscripts) : • ε12...n = 1 • for any index swap, ε changes sign: ε..i..j.. = - ε..j..i.. • therefore, if any two indices are the same, ε = 0 ε..i..i.. = - ε..i..i.. = 0 (B.2.8) Fact : If a is a permutation of z 0, then εa1a2... an = Parity( a) = (-1)Sa ≡ εa . (B.2.9) Proof : Since indices a i represent a permutation of z0 , it takes S a swaps to get from εa1a2... an to ε123...n by the definition of ε. In dense notation, one could say εa = (-1)Sa εz0 . Using (B.2.9) in (B.2.7) then gives these two classic det(M) expressions, Theorem 2: det(M) can be represented in these two ways: det(M) = Σ a1a2... an εa1a2... an M1a1M2a2 ...... Mnan (B.2.10a) det(M) = Σ a1a2... an εa1a2... an Ma11Ma22 ...... Mann . (B.2.10b) 57 In dense notation, one could write the above as ( a is now a vector), det(M) = Σ aεa Mz0a ( B . 2 . 1 0 a ) ' det(M) = Σ aεaMaz0 . (B.2.10b)' In most applications, the following simpler notation suffices: det(M) = Σ abc..q εabc..q M1aM2bM3c....Mnq ( B . 2 . 1 1 a ) det(M) = Σ abc..q εabc..q Ma1Mb2Mc3....Mqn . (B.2.11b) Theorem 3: For a square matrix M, adding a multiple of one row to another does not change det(M). The same is true for adding a multiple of one column to another. (B.2.12) Proof : (rows) Suppose we replace r3 → r3 + α r2. This says M 3i → M3i + α M2i . Eq (B.2.11a) says : det(M') = Σ abc..q εabc..q M1aM2b(M3c + α M2c) ....Mnq = det(M) + α Σ abc..q εabc..q M1aM2bM2c ....Mnq . Since M 1aM2bM2c ....Mnq is symmetric under b ↔c while εabc..q is anti-symmetric, the extra α term vanishes. That is to say, sum = Σbc AbcSbc = Σcb AcbScb = Σcb (-Abc)(Sbc) = - sum = 0. Using (B.2.11b) this argument shows that c3 → c3 + α c2 similarly does not alter det(M). Theorem 4: Swapping two rows (columns) of square matrix M causes det(M) → - det(M). (B.2.13) Proof : Let's swap rows 1 and 3 in M to get M'. Then from (B.2.11a), det(M') = Σ abc..q εabc..q M3aM2bM1c....Mnq = Σ abc..q [- εcba..q] M3aM2bM1c....Mnq // swap a↔c on ε = - Σ abc..q εabc..q M3cM2bM1a....Mnq // dummy rename a ↔c = - Σ abc..q εabc..q M1aM2bM3c....Mnq // reorder product = - det(M) . Using Theorem 1, we then know that det(M' T) = - det(MT), so swapping columns 1 and 3 negates det(M). An obvious corollary follows if the two sw apped columns have identical data : Corollary 4: If two rows (columns) of square matrix A are the same, then det(M) = 0. (B.2.14) 58 Here is one more well-known de terminant theorem which again demonstrates the power of the rearrangement theorem (B.2.4). Theorem 5 : det(AB) = det(A)det(B) (B.2.15) Proof : Since we already use A in a = Az0 , we shall prove det(XY) = det(X)det(Y) to avoid overloading symbols. Note that Parity( a) Parity( b) = (-1)Sa(-1)Sb = (-1)Sa+Sb = Parity(A b) from (B.1.4). Then : det(X)det(Y) = [ Σ a Parity( a) Π(z0,a; X)] [ Σb Parity(b) Π(z0,b; Y)] // (B.2.5) twice = Σa [ Σb Parity( a) Parity( b) Π(z0,a; X) Π(z0,b; Y) ] = Σ a [ Σb Parity(A b) Π(z0,a; X) Π(Az0,Ab; Y) ] // (B.1.6) with C = A = Σ a [ Σb Parity(A b) Π(z0,a; X) Π(a,Ab; Y) ] // A z0 = a = Σ a [ Σb Parity( b) Π(z0,a; X) Π(a,b; Y) ] // (B.2.4) with C = A = Σ b Parity( b) Σa Π(z0,a; X) Π(a,b; Y) // move Σa = Σ b Parity( b) Σa Π(z0,b; X Y ) = det(XY) . The idea here is to get a in the right place on both Π's so that Σ a Π(z0,a; X) Π(a,b; Y) = Σ a1a2... an X1a1X2a2 ...... Xnan Ya1b1Ya2b2 ...... Yanbn = Σ a1a2... an (X1a1Ya1b1)(X2a21Ya2b2).....(XnanYanbn) = (XY) 1b1 (XY)2b2 ... (XY) nbn = Π(z0,b; XY) which in dense notation one would write as Σ a Xz0aYab = (XY) z0b . Corollary 5 : det(ABC) = det(A)det(B)det(C) and so on. (B.2.16) Proof : det(ABC) = det(A[BC]) = det( A)det(BC) = det(A)det(B)det(C) . 59 B.3 Minors, Cofactors, and the Cofactor Expansions of det(M) Definition : The minor of a matrix element M rs is the determinant of the submatrix obtained by crossing out the rth row and the sth column. It is a challenge, however, to wr ite this out in symbols. Let's start with a more detailed version of (B.2.11a), where the a i are column summation indices, det(M) = Σ a1a2a3a4a5...an εa1a2a3a4a5...an M1a1M2a2M3a3M4a4M5a5 ...... Mnan . (B.3.1) We shall make the following conjecture for the form of minor(M 23) minor(M 32) = (-1)3-2 Σa1a2a4a5...an εa1a22a4a5...an M1a1M2a2M4a4M5a5 ...... Mnan (B.3.2) Compared to det(M) shown in (B.3.1), we have made these changes : • Removed the factor M 3a3 ( since row-3 matrix elements cannot appear in minor(M 32) ) • Removed the sum over a 3. • Replaced a 3 by the number 2 on the ε tensor. • added a sign factor (-1)3-2. Notice that since there is a 2 on the ε tensor, any time a summation index a i = 2 there is no contribution since then the ε tensor has two indices the same, so in effect the value 2 has been removed from all the residual summations a i. That is good, since column 2 is supposedly "crossed out" in minor(M 32). The factor (-1)3-2 is added so that the "diagonal term" in the minor will be positive. This sign (-1)3-2 gets used up if we slide the "2" to its natural position (position 2) on the ε tensor, minor(M 32) = Σa1a2a4a5...an εa12a2a4a5...an M1a1M2a2M4a4M5a5 ...... Mnan . (B.3.3) The diagonal term in minor(M 32) is then positive, matching the mechanical crossing-out method, ε 12345..n M11M22M44M55 ...... Mnn = + M11M22M44M55 ...... Mnn . (B.3.4) A more compact notation for (B.3.2) would be minor(M 32) = (-1)3-2 Σai,i≠3 εai,a3=2 Πi≠3 (Mi,ai) . (B.3.5) Starting over, we could show similarly that minor(M 42) = (-1)4-2 Σai,i≠4 εai,a4=2 Πi≠4 (Mi,ai) (B.3.6) We now have a sign factor (-1)4-2 because the "2" on ε has to be slid 2 positions to get to its natural location (position 4 on ε). More generally we can say minor(M s2) = (-1)s-2 Σai,i≠s εai,as=2 Πi≠s (Mi,ai) (B.3.7) 60 where the slide is now s-2 places. Still more generally we find, replacing 2 by r, minor(M sr) = (-1)s±r Σai,i≠s εai,as=r Πi≠s (Mi,ai) . (B.3.8) This then is our "best form" expression for a minor of matrix element M sr. Either sign will do, since (-1)s-r = (-1)s-r(-1)2r = (-1)s+r since (-1)2r = 1 . Choosing the + sign in (B.3.8) and putting (-1)s+r on the left side we get, (-1)s+r minor(M sr) = Σai,i≠s εai,as=r Πi≠s (Mi,ai) . (B.3.9) The left side here is called the cofactor of Msr, written cof(M sr). Thus we have shown that cof(M sr) ≡ (-1)s+r minor(M sr) = Σai,i≠s εai,as=r Πi≠s (Mi,ai) . (B.3.10) Fact : Neither minor(M sr) nor cof(M sr) are functions of the M sr matrix elements of M! In (B.3.10) this is so because: (1) row s is excluded in Π i≠s(Mi,ai); (2) ai = r is excluded by the factor ε ai,as=r . More intuitively, this is so because to get minor(M sr) we "cross out" row s and column r. Thus, ∂cof(Msr) ∂Msr = ∂minor(M sr) ∂Msr = 0 for any pair r,s in 1,2...n . (B.3.11) Finally, suppose we replace r with an integer which we call a s. Doing this causes εai,as=r → εai,as=as = εai = εa1a2a4a5...an = the normal ε tensor form . (B.3.12) Then (B.3.10) becomes the following, cof(M sas) ≡ (-1)s+as minor(M sas) = Σai,i≠s εai Πi≠s (Mi,ai) . (B.3.13) Now start again with det(M) of (B.3.1) and rewrite it in our compact notation: det(M) = Σ a1a2a3a4a5...an εa1a2a3a4a5...an M1a1M2a2M3a3M4a4M5a5 ...... Mnan = Σ ai εai Πi(Mi,ai) // next, extract the a s sum and its M sas factor : = Σas Msas [Σai,i≠s εai Πi≠s (Mi,ai)] // next use (B.3.12) to get = Σas Msas [cof(Msas)] // next, change summation index from a s to n = Σn Msn cof(Msn) . // valid for any s in 1,2...n (B.3.14) 61 This is the cofactor expansion of det(M) where one "works across row s" and n is a column index. Another form of this expansion is det(M) = det(M T) = Σn MT sn cof(MT sn) = Σ n Mns cof(Mns) ( B . 3 . 1 5 ) and in this form one "works down column s" wher e n is now a row index. We have just proven : Theorem 6 (Cofactor Expansions): det(M) = Σ n Msn cof(Msn) = Σn (rs)n cof(Msn) work across row s s = 1,2....n (B.3.16a) det(M) = Σn Mns cof(Mns) = Σn (cs)n cof(Mns) work down column s s = 1,2....n (B.3.16b) The term "cofactor" presumably arose because each sum term consists of a "factor" M ij and a "co-factor" which cooperates with the factor to make the su m, as coauthors cooperate to write a book. A common notation is to define Mij ≡ cof(Mij) and then the cofactor expansions become det(M) = Σ n MsnMsn det(M) = Σn MnsMns . Although compact, this is not so handy for tensor notation with its up and down indices (see Lucht Ref). B.4 Expressions for the inverse matrix M-1 and Cramer's Rule Consider the expansion (B.3.16a) working across row s, det(M) = Σ n Msn cof(Msn) . ( B . 4 . 1 ) Suppose we replace the elements of row s (M sn) with the elements of some other row r in M (M rn). In doing so we have created a new matrix, call it M'. Notice that cof(M' sn) = cof(M sn) because , although going from M to M' we have altere d row s, we have not altered cof(M sn) since this depends only on the entries in all the rows other than row s (think "crossing out" row s for minor(M sn) ). See Fact (B.3.11) and text above. Therefore we find, det(M') = Σn M'sn cof(M'sn) = Σn Mrn cof(Msn) . (B.4.2) But by Corollary 4 (B.2.14) det(M') = 0 si nce two rows of M' are the same. Thus 0 = Σn Mrn cof(Msn) . r ≠ s ( B . 4 . 3 ) Combining this with (B.4.1) gives, Σ n Mrn cof(Msn) = det(M) δr,s . (B.4.4) 62 Now for clarity define a cofactor matrix C in this manner C sn ≡ cof(Msn) . ( B . 4 . 5 ) In matrix notation, we could define the matrix cof(M) to be matrix C and then C = cof(M) ⇒ Csn = [cof(M)] sn = cof(M sn) . (B.4.6) Now (B.4.4) says Σ nMrn Csn = det(M)δ r,s or Σ nMrn CT ns = det(M) δr,s and finally in matrix notation, MCT = det(M) 1 or M [ CT det(M) ] = 1 . (B.4.7) Therefore we find this classic square matrix inversion formula, Theorem 7: M -1 = CT det(M) = [cof(M)]T det(M) = cof(MT) det(M) (B.4.8) To verify the last equality in (B.4.8) consider : = [cof(M)] T ns = [cof(M)] sn // meaning of transpose = cof(M sn) // (B.4.6) = cof(M T ns) = [cof(MT)]ns // (B.4.6) applied to MT and therefore we have this matrix identity, [cof(M)] T = [cof(MT) ] . ( B . 4 . 9 ) Using (B.4.8) it is trivial to solve a non-singular (detM ≠0) system of linear equations y = M x : y = M x ⇒ x = M-1 y = 1 det(M) cof(MT) y = 1 det(M) [cof(M)]T y . (B.4.10) 63 In components, x s = 1 det(M) Σn ([cof(M)]T)sn yn = 1 det(M) Σn yn [cof(M)] ns = 1 det(M) Σn yn cof(Mns) . (B.4.11) Now recall the cofactor expansion (B.3.16b), det(M) = Σ n Mns cof(Mns) = Σn (cs)n cof(Mns) . (B.4.12) If one replaces column c s in M by y one gets det(M[ c s→ y]) = Σn yn cof(Mns) // (B.4.12) = det(M) x s . // (B.4.11) (B.4.13) Solving this for x s we obtain another classic result, Theorem 8 (Cramer's Rule, 1750): // Gabriel Cramer (1704-1752), Swiss, made it to age 47 y = M x ⇒ x s = det(M[ cs→ y]) det(M) ( B . 4 . 1 4 ) Comment: Names of Matrices Historically the matrix CT = [cof(M)]T = [cof(MT)] was called the adjoint of matrix M and was often denoted by M^ and then (B.4.8) reads M-1 = M^ /|M|. But this then conflicted with another usage, namely the matrix M† ≡ (MT)* = (M*)T being the adjoint of M (the conjugate transpose of M). When M† = M, M is said to be self-adjoint and has significance in quantum (matrix) mechanics: the quantum operators of all physical observables are self-adjoint and thus have real eigenvalues. Nowadays the matrix [cof(M)] T is called the classical adjoint of M, while M† is then the Hermitian adjoint (or Hermitian conjugate ) of M and is sometimes written MH. If M† = M then M is said to be Hermitian. Historically when [cof(M)]T was the adjoint of M, M† was called the associate of M. The classical adjoint sometimes appears as the adjugate or adj(M). Transpose matrices MT are often denoted M~ . On a related matter, some older authors used the word minor to refer to what we now call a cofactor. See for example Morse & Feshbach page 509. These authors never use the word cofactor. 64 References The method of Lagrange multipliers often makes cam eo appearances in textbooks on the calculus of multiple variables. There are also several pages and pdf's on the web, search on "Lagrange Multipliers". Whole books on the use the Lagrange multipliers for sp ecialized applications may be found at the usual book vendor websites. R.C. Buck with E.F. Buck, Advanced Calculus , 2nd Ed. (McGraw-Hill, New York, 1965). The discussion of extremum problems with constraints occupies pp 359-363. In the 3rd Edition (McGraw-Hill, New York, 1978), reprinted by (Waveland Press, Long Grove IL, 2003), this discussion has moved to pages 536-540. This excellent and now classic book was first published in 1956. H. Goldstein, Classical Mechanics (Addison-Wesley, Boston, 1950), another classic. There is a 3rd edition 2001 by Goldstein, Safco and Poole. S. Jensen, An Introduction to Lagrange Multipliers, http://www.slimy.com/~steuard/teaching/tutorials/Lagrange.html ( includes video lectures) P. Lucht, Tensor Analysis and Curvilinear Coordinates (2014), http://user.xmission.com/~rimrock . If not there, search on . P.M. Morse and H. Feshbach, Methods of Theoretical Physics ( McGraw-Hill, New York, 1953). G.E. Shilov, Linear Algebra (Dover Publications, 1977). W. F. Trench, The Method of Lagrange Multipliers (2013), http://digitalcommons.trinity.edu/mono/7/ , Supplement 2 (31p). This paper also has a "Theorem 1" which is related to but different from ours. wiki on Lagrange Multipliers: https://en.wikipedia.org/wiki/Lagrange_multiplier