MATH/APPM 4650
Class notes, May 1
Chapter 6 review

This chapter covers methods for solving the matrix equation A x = b for x, where A is an n × n square matrix and b is a given n-dimensional column vector.

This sort of problem comes up frequently in many applications, including curve-fitting (splines) and regression, boundary-value problems for ODEs, implicit methods for PDEs, and solving systems of nonlinear equations.

Our goal is to find algorithms that work well as n gets very large. We also want to minimize roundoff error.

• For a square matrix, the system has a unique solution if A is nonsingular. Equivalent conditions for this are that the rank of A is n, that A can be row-reduced to an upper triangular matrix with all diagonal entries nonzero, that the determinant of A is not zero, and that there is an inverse matrix A^-1.
• (5.1) The fastest general-purpose algorithm is Gaussian elimination, which reduces the matrix to upper triangular form. It involves three nested loops (loop on the columns you want to clear, loop on the entries of each column, loop on the remainder of each row so you don't change the equation). This algorithm takes n³/3 + n²/2 - 5n/6 multiplications/divisions and n³/3-n/3 additions/subtractions. The important thing to remember is the highest term: roughly n³/3 operations no matter how you count them.
• After reducing to upper triangular form, the remainder of the algorithm is backward substitution. It requires roughly n²/2 operations, since it's two nested loops (one loop for each component of x, the second loop since successive x values are computed as a sum of previous values).
• Gauss-Jordan elimination, which reduces the matrix to diagonal form, requires roughly n³/2 operations. It is definitely slower in every case than Gaussian elimination.
• (5.2) The basic Gaussian elimination algorithm does not take into account the possibility that some diagonal element will be zero or small. Since it requires dividing by all diagonal elements, one should incorporate a pivoting algorithm to switch rows if a diagonal element is too small. Properly done, such an algorithm takes roughly O(n²) operations which is insignificant compared to n³/3.
• We can either do partial pivoting (pick out the highest absolute value of a term in any column) or scaled partial pivoting (pick out the highest absolute value of the ratio of a row to its largest element). Notice that the scale factors must be computed at the beginning of the procedure, not line by line, or else you add an O(n³) operation to the reduction. Scaled partial pivoting is equivalent to multiplying every equation by a factor at the beginning so that the largest coefficient in every equation is one, then doing ordinary partial pivoting on what's left.
• (5.3) The inverse of a matrix, if it exists, may be found by applying Gaussian elimination, with vector b replaced by the identity matrix. It requires roughly 4n³/3 operations this way. Factorization is usually faster.
• (5.4) Determinants are generally very difficult to compute directly; expanding on rows or columns without taking advantage of special features will take about n! operations, which grows much faster than even an exponential aⁿ. For large matrices it is much better to do Gaussian elimination first: Gaussian elimination never changes the determinant, and Gaussian elimination with partial pivoting will multiply the determinant by +1 or -1 depending on how many rows you switch. This takes roughly n³/3 operations, just like Gaussian elimination itself.
• (5.5) The best way to solve the equation Ax = b repeatedly is to first find a matrix factorization A = LU, where L is lower triangular and U is upper triangular. Then you solve LU x = b by first solving L y = b for the auxiliary vector y (which takes O(n²) operations) and then solving U x = y for the desired vector x (which also takes O(n²) operations).
• The standard way to get this is to assume L has ones on its diagonals. Then U is the same matrix you'd get from Gaussian elimination, and L is the matrix whose off-diagonal elements are the multipliers for Gaussian elimination.
• If you need to switch rows, you can incorporate this via a permutation matrix. A = P^tLU.
• (6.6) If the matrix is symmetric, then the factorization is A = LDL^t for a diagonal matrix D. If you already have LU, just pull out the diagonal elements from U and divide the matrix through.
• If the matrix is strictly diagonally dominant (the sum of the abs. vals. of the off-diagonal elements in each row is strictly less than the abs. val. of the diagonal element), then the matrix is guaranteed to not be singular and also to never need row switching in Gaussian elimination.
• If the matrix is positive definite (symmetric and all eigenvalues are strictly positive), then D in the LDL^t factorization is positive, and its square root exists, so you can write A = L'L'^t where L' = Lsqrt(D).
• If the matrix is a band matrix (e.g. tridiagonal) then Gaussian elimination takes only O(n) operations, since you don't need to extend the loops all the way out. This is a huge time savings.