Jacobians of Matrix Transforms

Chapter 8 Jacobians of Matrix Transforms (without wedge products) Chapters 8 and 9 provide the background to derive the matrix Jacobians such as those found in Chapter 42. 8.1 Matrix and Vector Differentiation In this section, we concern ourselves with the differentiation of matrices. We begin with the familiar product rule for scalars, d(uv) = u(dv)+v(du), from which we can derive that d(x)=3zdr. We refer to dz as a differential. We all unconsciously interpret the "dz" symbolically as well as numerically. Sometimes it is nice to confirm on a computer that (+)-32 (8.1) I do this by taking to be 1 or 2 or randn(1) and to be .001 or 0001. The product rule holds for matrices as well: d(UV)U(dV) + (AU)V. In the examples we will see some symbolic and numerical interpretations. Example 1: Y = X We use the product rule to differentiate Y(X) X³ to obtain that d(X)=X2(dX) +X(AX)X+(AX)X2. When I introduce undergraduate students to matrix multiplication, I tell them that matrices are like scalars, except that they do not commute. The numerical (or first order perturbation theory) interpretation applies, but it may seem less familiar at first. Numerically take X-randn(n) and E-randn(n) fore.001 say, and then compute (X+E)-X XE + XEX + EX2. (8.2) This is the matrix version of (8.1). Holding X fixed and allowing E to vary, the right-hand side is a linear function of E. There is no simpler form possible. Symbolically (or numerically) one can take dXE which is the matrix that has a one in element (k,I) and 0 elsewhere. Then we can write down the matrix of partial derivatives: =X2(Ext)+X(E)X+(E)X2. In general, the directional derivative of Yiy(X) in the direction dX is given by (dY). For a particular matrix X, dY(X) is a matrix of directional derivatives corresponding to a first order perturbation in the direction EdX. It is a matrix of linear functions corresponding to the linearization of Y(X) about X. Structured Perturbations We sometimes restrict our E to be a structured perturbation. For example if X is triangular, symmetric, antisymmetric, or even sparse then often we wish to restrict E so that the pattern is maintained in the perturbed matrix as well. An important case occurs when X is orthogonal. We will see in an example below that we will want to restrict E so that. XTE is antisymmetric when X is orthogonal. Example 2: y Here y is a scalar and dot products commute so that dy2rds. When y 1, x is on the unit sphere. To stay on the sphere, we need dy = 0 so that rdr=0, i.e., the tangent to the sphere is perpendicular to the sphere. Example 3: y=Ar Again y is scalar. We have dy dr Ar+z7 Adz. If A is symmetric then dy = 227 Adr. Example 4: Y=X-1, We have that XYI so that X (dY) + (dX)Y=0 so that dYXaxx. We recommend that the reader compute e ((X+E)-1-X-1) numerically and verify that it is equal to -X-1EX-1, Example 5: 1 QTQ. If Q is orthogonal we have that QdQ+dQQ=0 so that QdQ is antisymmetric. Notation: Let () T ER. We denote a vector of differentials by dr = dan dr 8.4 Jacobians of Matrix Factorizations (without wedge products) Let A be an nxn matrix. In elementary linear algebra, we learn about Gaussian elimination, Gram-Schmidt orthogonalization, and the eigendecomposition. All of these ideas may be written compactly as matrix factorizations, which in turn may be thought of as a change of variables: Here is a table that is expanded in Section 35.1. Gaussian Elimination: A= L- ↑ unit lower triangular parameter count U ↑ upper triangular Gram-Schmidt: A= Q ↑ R ↑ Orthogonal upper triangular n(n-1)/2 + n(n+1)/2 n(n-1)/2+n(n+1)/2 Eigenvalue Decomposition: A= X ↑ AX-1 (nn) + n eigenvectors eigenvalues ↑↑ eigenvector eigenvalue Each of these factorizations is a change of variables. Somehow the n² parameters in A are transformed into n² other parameters, though it may not be immediately obvious what these parameters are (the n(n-1)/2 parameters in Q for example). Our goal is to derive the Jacobians for those matrix factorizations. 8.4.1 Jacobian of Gaussian Elimination (A = LU) In numerical linear algebra texts it is often observed that the memory used to store A on a computer may be overwritten with the n² variables in Land U. Graphically the nxn matrix A U A L Indeed the same is true for other factorizations. This is not just a happy coincidence, but deeply related to the fact that n² parameters ought not need more storage. Theorem 8.1. If ALU, the Jacobian of the change of variables is det J.-10-1 = П (=1 Proof 1: Let ALU, then using dA = LAU + ALU L-HAU- = QUU+LAL (UL)dA (UTI)-dU+ (IL)-dL or dA = (UL) ((UI)-dU+ (IL)-dL). -2 Using the triangular case of Section 8.3, the Jacobian of (U) is ... and that of UTL is (det U) un putting the two together we obtain that det That is the fancy slick proof. Here is the direct proof. It is of value to see both. Proof 2: Let Mij as in the diagram above. Since ALU, we have that MikM+ MkM+ for ij. A1 if isj Misty if ij for i A and Aj = Therefore (8.6) Notice that Aij never depends on Mag when pi orqj. Therefore if we order the variables first by row and next by column, we see that the Jacobian matrix is lower triangular with diagonal entries given in (8.6). Remembering that the determinant of a lower triangular matrix is the product of the diagonal entries, the theorem follows. Π Perhaps it should not have been surprising that the condition that A admits a unique LU decomposition may be thought of in terms of the Jacobian. Given a matrix A, let pa denote the determinant of the upper left k by k principal minor. It may be readily verified that ut=P and hence the Jacobian is pip2P-1. The condition that A admits a unique LU decomposition is well known to be that all the upper-left principal minors of dimension smaller than a are non- singular. Otherwise the matrix may have no LU decomposition. For example, (1) has no LU factorization. It may also have many LU decompositions as does the zero matrix when > 1. This degeneracy ex- plains the need for pivoting strategies (ie., strategies that reorder the matrix rows and/or columns) for solving linear systems of equations even if the computation is done in exact arithmetic. Modern. Gaussian elimination software for solving linear systems of equations in finite precision include pivoting strategies designed to avoid being near a matrix with such a degeneracy. Exercise. If A is symmetric positive definite, it may be factored ALLT. This is the famous Cholesky decomposition of A. How many independent variables are in Af In L? Prove that the Jacobian of the change of variables is det J 2 +1-4 Research Question. It seems that the existence of a finite algorithm for a matriz factorization is linked to a triangular Jacobian matriz. After all, the latter implies only one new variable need be substituted at a time. This is the essential idea of Doolittle's or Crout's matriz factorization achemes for LU.