
Software: 
Matrix Addition or Subtraction To find [C]=[A] ± [B] we obtain each element C_{ij} of [C] by C_{ij} = A_{ij} ± B_{ij} Matrix Multiplication To find [C]=[A][B] we obtain each element C_{ij} of [C] by
n where n = # of columns of [A] = # of rows of [B]. Determinant To find the determinant of matrix [A], we triangularize it using elementary row operations. The determinant is then the product of the diagonal elements. Trace The trace of matrix [A] is the sum of its main diagonal elements. Inverse The inverse of [A] is computed using the GaussJordan method. It starts by augmenting [A] on the right by an identity matrix [I] of the same order. Then, by a series of elementary row operations, we transform [A] into an identity matrix, making sure to extend each operation to include the entire augmented matrix. When the process is finished, the original [I] with which we flanked [A] will have been transformed into the inverse of [A]. Adjoint The adjoint of a matrix is the transpose of the matrix of the cofactors of its elements. If [B] is the adjoint of [A], then the elements of [B] are B_{ji} = (1)^{i+j} det [Mij] where [Mij] is a submatrix of [A] obtained by deleting from [A] the ith row and the jth column. If [A] is not singular, its adjoint may be found using the relation adj [A] = det [A] * [A]^{1} where adj [A] is the adjoint of [A], det [A] is the determinant of [A], and [A]^{1} is the inverse of [A]. Transpose To find the transpose of a matrix [A] we form a matrix [A]^{T} whose nth row is equal to the nth column of [A]. LU Factors To find the [L][U] factors of matrix [A], essentially what we do is carry out the process of Gaussian Elimination, in which a system of linear equations [A][X]=[B] is solved by uppertriangularizing [A] by the process of elimination, to get [U][X]=[C] and solving this system by backsubstitution. The [U] that appears in the last system above is the [U] factor we seek, while the elements of the [L] factor are the factors Fij by which row j is multiplied before subtracting it from row i during the elimination process. QR Factors The GramSchmidt orthogonalization process is used to find the QR factors. If we have a set of vectors {A_{1 }, A_{2 },...A_{n}} we can generate an orthonormal set {Q_{1},Q_{2},...Qn} as follows: Choose the first vector in the direction of A1:
V_{1} then continue thus:
V_{2}
V_{3} Now, since we want [R] such that [Q][R]=[A], we premultiply both sides by [Q]^{T}, which is the transpose of [Q]: [Q]^{T}[Q][R]=[Q]^{T}[A] but, since [Q] is orthonormal, then [Q]^{T}[Q] = [I], so [R]=[Q]^{T}[A] Test for Definiteness For illustration purposes, consider the 4th order symmetric matrix:
é
A_{11} A_{12} A_{13} A_{14}
ù where A_{ij}=A_{ji}. [A] is positive or negative definite if the product [x
y z t] é
A_{11} A_{12} A_{13} A_{14}
ù
éxù is positive or negative, respectively, for all nontrivial [x y z t]. Clearly, we can't test that product for all possible [x y z t], but there are other equivalent tests. One is to check the signs of the eigenvalues: if they are all positive or all negative, the matrix is positive or negative definite, respectively. Another test is to compute the determinants A_{1 }, A_{2 }, A_{3} and A_{4} of these principal submatrices: A_{1} = det [ A_{11} ] A_{2}
= det é
A_{11} A_{12} ù A_{3}
= det é
A_{11} A_{12} A_{13}
ù A_{4}
= det é
A_{11} A_{12} A_{13} A_{14}
ù [A] is positive definite if A_{i} > 0 for i = 1, 2, 3, 4 (all i) [A]
is negative definite if
A_{i} < 0
for
i = 1, 3 (i odd) The matrix may also be positive semidefinite, negative semidefinite, or indefinite. Simultaneous Linear Equations The method used is called Gaussian Elimination. Suppose the system is [A][X]=[B] where [A] is nbyn, [X] is nby1 and [B] is nby1. Here [A] is the matrix of coefficients, [X] the solution sought and [B] the righthandside of the system. We start by augmenting matrix [A] by attaching to it the column [B]. Then, by a series of elementary row operations, we reduce the extended matrix to echelon form, that is, [A] is made upper triangular, making sure to extend each operation to include column [B]. If we call [T] the triangularized [A] and [B'] the transformed [B], then we simply solve the system [T][X]=[B'] by back substitution. Simultaneous Linear Equations (Solution by [L][U] Factors) This option solves the equations using the LU factorization of the matrix of coefficients. This method is equivalent to the method of Gaussian Elimination, which actually consists first of elimination and then backsubstitution, but it facilitates the solution of many systems sharing the same matrix of coefficients by avoiding the duplication of many operations. Once the LU factors are known, the solution is reduced to a forwardsubstitution and then a backsubstitution. Suppose the system is [A][X]=[B] where [A] is nbyn, [X] is nby1 and [B] is nby1. Here [A] is the matrix of coefficients, [X] the solution sought and [B] the righthandside of the system. We start by factoring [A] into its LU factors. The above equation becomes [L][U][X]=[B] Each time we are given a set of independent terms [B], we first solve for [C] the equation [L][C]=[B] which, since [L] is lowertriangular, is easily solved by forwardsubstitution, then we solve for [X] the equation [U][X]=[C] which, since [U] is uppertriangular, is easily solved by backsubstitution. Of course, the factorization [A]=[L][U] may really be [A]=[P][L][U], where [P] is a permutation matrix, but this may be remedied by reorganizing the independent terms using the map of row exchanges. Overdetermined or Inconsistent Systems (LeastSquares Solution) If we have a system with more equations than unknowns, [A][X]=[B] where [A] is the matrix of coefficients, [X] the solution sought, and [B] the system's righthandside, we may not be able to find a solution which satisfies all the equations, but we can find the leastsquares solution by solving [A]^{T}[A][X']=[A]^{T}[B] where [A]^{T} is the transpose of [A] and [X'] is the leastsquares solution. Eigenvalues and Eigenvectors In principle, to obtain the eigenvalues of a matrix its characteristic equation may be solved: det [A]  l[I] = 0 where [I] is the identity matrix. For an nbyn [A], it is of nth degree in l. Although numerically extracting the roots of an nth degree equation may be relatively easy, the expansion of the determinant is impractical even for a 4by4 matrix. For large matrices it may be prohibitive, hence the need for algorithms to extract eigenvalues. The shifted QR algorithm generally converges relatively fast. If we could make [A] uppertriangular while preserving its eigenvalues, we would find those eigenvalues in its main diagonal. Elementary row operations destroy the eigenvalues, so they can't be used, and generally [A] can't be triangularized in a finite number of steps. What is done is to make [A] as nearly upper triangular as possible while preserving its eigenvalues. This is accomplished by the QR method, except for 2by2 blocks which correspond to complex conjugate pairs of eigenvalues. Basically, at each step we factor [A] into [Q][R] using the GramSchmidt orthogonalization process. The factors are then reversed to obtain a more nearly triangular [A'] with the same eigenvalues. In practice, the shifted QR algorithm is the one used, which provides faster convergence. Furthermore, [A] is first transformed to produce as many zeroes as possible. It is made upperHessenberg, which is nearly uppertriangular except for the elements of the diagonal below its main diagonal. This is done using Householder transformations. This preliminary step insures faster QR factorization. It should be said that the upperHessenberg form is preserved at each step of the QR process. To illustrate, call [A_{0}] the original matrix already made upperHessenberg. Then the QR factorization will give us [A_{0}] = [Q_{0}][R_{0}] The next [A] is, reversing the factors: [A_{1}] = [R_{0}][Q_{0}] Again we factor [A_{1}], then reverse the factors: [A_{1}] = [Q_{1}][R_{1}] ; [A_{2}] = [R_{1}][Q_{1}] and so on. However, in the shifted algorithm, if l_{k} is close to an eigenvalue at step k, we shift [A_{k}], so [A_{k}]  l_{k}[I] = [Q_{k}][R_{k}] where [I] is the identity matrix. At each step [A] is made more and more upper triangular as the elements below its main diagonal approach zero (except, as was said before, for those which correspond to complex eigenvalues). The process is stopped when all the eigenvalues have converged, or when it has gone through the number of cycles specified. An eigenvalue is assumed to have converged when its value differs from its previous value by less than the specified error. In practice, the first one to converge is frequently the smallest one. If this is the one sought, a moderate number of cycles may produce the desired result. The eigenvectors can be found using the Shifted Inverse Power algorithm. It starts by making [A] as nearly singular as possible by shifting it using the eigenvalue for which the eigenvector is sought. The iteration for l_{k} starts with: [ [A]  l_{k}[I] ][X] = [RND] where [RND] is a column vector with random components. The system is solved and the solution is then used in place of [RND] to obtain the next better [X]. It cycles until the vector solution changes very little, that is, until each component of the normalized eigenvector changes by less than the specified error. The above computation of eigenvectors allows for the possibility of the existence of several distinct eigenvectors corresponding to an eigenvalue which is a multiple root of the characteristic equation. If the successive approximations sequence were to be always started with [1] or any other fixed column vector instead of [RND] then it would always converge to the same eigenvector. 

Copyright © 20012010 Numerical Mathematics. All rights reserved.
