Lecture 10: Linear Systems (Matrix Magic)

Matrices and Linear Systems in Attaway

The following lecture overheads and accompanying reading are independent of MatLab. Attaway Chapter 11.1 links matrices as described here with their representation and use in Matlab. Attaway 11.1 should be read and understood as a complement to the on-line lecture and reading.

Attaway 11.2, leading up to 11.2.2.1, gives the background and algorithm for Gaussian Elimination as a method of solving Linear Equations. It is useful "how-to" reading for the assignment.

Linear Systems

• A linear equation is one which contains only scalar multiples of its variables and constants. For example:
  a₁X₁ + a₂X₂ + a₃X₃ = c
  is a linear equation with three variables.
  X₁-X₃ represent the variables, a₁-a₃ are the corresponding coefficients, and c is called the constant term.
  
  5x + 7.5y - 10z = 30
  
  is a particular example.
• A system of linear equations is a (finite) set of such equations that are (generally) linked by containing (some of) the same variables.
• Linear systems are important because they are useful for modeling a large number of practical situations and mathematicians know a (heck-of-a) lot about them, including:
  • Exactly when solutions exist.
  • Exactly when a solution is unique.
  • Efficient algorithms for finding solutions when they exist.
  • Effective methods for approximating solutions in cases where actual solutions do not exist (least squares methods).
• In general, if we have N variables (or unknowns), we need N linearly independent equations to find a unique solution. Such independence means we cannot derive any of the equations from the others.
• If we have more variables than equations, the system is said to be underdetermined. The equations will generally constrain the solution to a linear subspace of the space of possible solutions, but there is no single, unique solution. E.g. with three variables, a single linear equation describes a plane (subspace of 3-D) each of whose points satisfies the equation.
• If we have more equations than variables, then the system will, in general, have no solution (unless some of the equations are linearly dependent). Such a system is said to be overdetermined or inconsistent. Although there may be no actual solution, there are often points in the space of variables that are "almost solutions" in a sense that can be made mathematically rigorous. Least squares analysis is one common approach to finding such approximate solutions. (More in the Model Fitting topic coming up).

Geometric Intuition

• Recall that the points satisfying the linear equation
  ax + by = c
  in the two variables x and y, fall on a line in the x-y plane.
• If we introduce a second equation with (generally) different constants a, b, and c, then we have a second line. There are four possibilies:
1. In general, a second line will intersect the first in a single point. This is the generic case of two unknowns, two equations, and a unique solution.
2. If the two lines are parallel, (e.g. with a and b the same and c different), then there is no solution. This is an example of an inconsistent system.
(Also if there are three or more lines that do not intersect at a single point. This case is often referred to as overdetermined)
3. If the two lines are identical, (e.g. with a, b, and c in the second equation all the same multiple of their values in the first equation) then there are an infinite number of solutions (all the points on the line). In this case, the equations are linearly dependent in a particularly simple way. This is an example of an underdetermined system.
4. If a, b, and c are all 0 in both equations, then the system is said to be trivial, and all points are solutions. This is can be viewed as an extreme example of an underdetermined system.
• In three dimensions, equations represent planes, and analogous geometric intuitions hold (e.g. three planes generally intersect in a point. The opportunities for underdetermined systems are more complex than coinciding planes (e.g. three planes can intersect in a line). The opportunities for inconsistent systems are also more complex than parallel planes (e.g. three planes can intersect pairwise in three parallel lines).

Matrix Representation

• A system of linear equations can be represented compactly using a matrix of the coefficients and (column) vectors for the variables and constant terms.
  Quick intro to matrices
• The multiplication of a vector by a matrix is defined so that performing the manipulation generates the linear equations.
• For example, suppose we have the following system of equations:
  a₁x + b₁y = c₁
  a₂x + b₂y = c₂
  This can be written in matrix form as:

[a1 b1] [x] = [c1]
[a2 b2] [y]   [c2]

• More generally, we will use single variable names and subscripts to refer to matrix and vector elements, where the elements of matrix A are referred to as a_{row col}.
  The following thus represents a generic 3x3 system.

[a11 a12 a13] [x1]   [c1]
[a21 a22 a23] [x2] = [c2]
[a31 a32 a33] [x3]   [c3]

• Compactly, in matrix-vector notation, the above is written:
  Ax = c
• There are Elementary Row Operations (ERO's) that produce a "new" system of equations that has the same set of solutions as the original and may be more fit for computation:
  • Swapping ERO: the rows of the matrix can be interchanged or rearranged without changing the system of equations represented as long as the elements of the constant c vector are rearranged in the same way.
  • Linear Combination ERO: a multiple of a row (and its corresponding element in the constant vector) may be multiplied by a constant, or added to another row, or both, without changing the underlying system (its solutions).
• Similarly, the columns of the matrix can be rearranged without changing the system of equations represented, as long as the elements of the variable x vector are rearranged in the same way.
• Systems with more or fewer equations than variables can be written in (non-square) matrix-vector form. In this case the number of columns of the matrix equals the number of variables (and the length of the variable vector), and the number of rows of the matrix equals the number of equations (and the length of the constant vector).

Linear Systems Happen

• Mechanical Force Analysis
• LCR Circuit Analysis
• Optics: (e.g. Formalize Ray Tracing): 160 Ex.
• Chemistry and Chem. Engg: Reaction Kinetics
• General: Fitting Mathematical Models to Data: 160 Ex.
• General: Least Squares Optimization
• General: Numerical Solution of Differential Equations: 160 Ex.

Gaussian Elimination

• Gaussian elimination is a general algorithm for solving systems of linear equations
• Seems to have been known by the Ancient Chinese prior to 100 BCE.
• European introduction by Carl Friedrich Gauss (German Mathematician) around 1809 in context of least squares analysis.
• Basic technique, with minor variations, is still used for general systems up to several thousand variables.
• For large systems (i.e millions of variables) iterative approximation techniques are used.
• For sparse systems (mostly 0s in the matrix) a variety of special, fast methods have been developed.

Basic Approach

• Recall the linear combination ERO: Adding a multiple of one equation to another produces a new equation that holds iff the original equations hold.
• Idea is to add a multiple of one equation to another so that the coefficient of some (selected) variable is 0.
• The result has one less variable than the original.
• By doing this progressively on selected variables, in an organized fashion, systems of equations with fewer and fewer variables can be produced, until an equation with only one variable is obtained.
• This last equation can be solved directly for the last variable, which can be substituted into the 2-equation system to obtain the value of a second variable, and so forth until values for all the variables have been obtained.

Details A: Reduction Step

Implementation Note: Though it mixes element semantics and is not mathematically a pure object, the augmented matrix of the N-variable system AX = B is the
N x N+1 matrix [A | B] (that is the column of constants B concatenated onto the right of A ). This augmented matrix is useful because all the EROs involve both A and B , which are now together in the augmented matrix and can be manipulated together twice as easily as apart.

• Equations are written in Ax = c matrix form. The x vector is sometimes not written out at every step as it serves only to specify the order of the unknowns.
• Appropriate multiples of the first row are added to the other rows so that the first coefficient is 0 in each subsequent row. This produces a column of 0s below the (1,1) element of the matrix. The constant elements are treated as part of the row.
• The appropriate multiples are determined by dividing the first coefficient of each lower row by the first coefficient of the first row (the (1,1) element). This divisor is known as the pivot
• In similar fashion, appropriate multiples of the second row are added to the rows below it to produce 0s in the second column below the (2,2) element (which now serves as the pivot)
• The process is repeated for subsequent rows until an upper triangular matrix that contains 0s below the main diagonal is produced. This matrix, along with the (modified) constant vector, represents a system that has exactly the same solutions as the initial system. This completes the reduction stage.

Details B: Back Substitution Step

• The algorithm now enters the back substitution stage. Note that the last row of the upper triangular matrix represents an equation in one variable (the one associated with the last column), which can be directly solved for that variable.
• The next-to-last row represents an equation in two variables, one of which is the variable just solved for. Substituting in the value, yields an equation that can be solved directly for the variable associated with the next-to-last column.
• The process is repeated until values for all the unknown variables have been obtained. This completes solution of the system.
• Worked example

Computational Complexity

• If we count up the operations involved, it turns out that approximately N³/3 additions and multiplications are needed for the reduction step, and approximately N²/2 additions and multiplications for the back-substitution step, where N is the size of the system.
• As N becomes large, the operation count is dominated by the N³/3 term. Computer scientists describe the situation by saying that the algorithm is order of N³.
• This is conventionally written O(N³), referred to as Big O notation. The technical meaning is that for sufficiently large N the operation count is bounded by kN³ for some constant k.
• Such asymptotic bounds are an important area of study in the field of computer science known as computational complexity theory.
• Rather surprisingly, kN³ is not a lower bound. Complexity is the same as square matrix multiplication for which (mostly impractical) algorithms of lower asymptotic complexity are known, e.g., Strassen's Algorithm which is approximately O(N^2.807).

Potential Problem A: Zero Pivot

• If the diagonal element that is to serve as a pivot is 0, no multiple of the row can eliminate that variable in the following rows.
• An easy solution is to use the swapping ERO to swap the the problematic row (and the associated constant term) with a lower row that does not have a 0 in the column being reduced.
• If no such row exists, then in one sense, the column is already reduced and we can go on to the next one. However we have discovered a linear dependency, which means that a unique solution does not exist. Such a system is said to be singular.
• More on reduction in singular systems

Potential Problem B: Near Zero Pivot

• If a pivot element is very small compared to one or more of the elements below it in the column, then a large multiple of the pivot row must be added to bring the column coefficient to zero. This will amplify any errors existing due to roundoff or measurement error, and can cause an inaccurate solution of the system.
• A simple approach is the swap rows to use the one that currently has the largest value in the pivot column. This is known as partial pivoting.
• A more complex approach is to swap both rows and columns so as to obtain the largest possible pivot. This is known as full pivoting, but is not generally used as the search is time consuming, and partial pivoting is usually sufficient.

Potential Problem C: Unstable System

• If a situation arises where the best available pivot is small compared to values above it (in the rows where reduction has already been completed) then large multiplications will be introduced during the back-substitution stage, which can also amplify existing errors.
• In this case, the system is intrinsically sensitive to error. Such a system is termed ill-conditioned or sometimes nearly singular.
• An example is a pair of equations in the plane representing nearly parallel lines. Clearly a small change in the position or orientation of one of the lines can cause a large change in the location of the intersection point (or even make them parallel).
• An ill-conditioned system arising from real experimental data often indicates a deeper problem with experimental design, data acquisition protocols, or even with the underlying model. Detecting such systems is thus important.

Variations

• Instead of proceding with the back-substitution step after obtaining an upper triangular form, the reduction formality can be continued, using the last row to eliminate the last column coefficients, the next-to-last row to eliminate the next-to-last column coefficients, and so forth. This produces a diagonal form that represents a separate, easily solved equation for each variable. This is known as Gauss-Jordan elimination.
• Each row and corresponding constant of the diagonal form can be divided by the value of the diagonal element, producing the identity matrix, from which the values of the variables can be simply read off.
• If we place the identity matrix I adjacent to the original coefficient matrix A instead of the constant vector c, and carry out the above reduction of A to the identity, adding multiples of the rows all the way across, the identity matrix is transformed into A^-1. This is a standard algorithm for computing the inverse.