MATH/APPM 4650
Class notes, April 13
Section 6.1: Linear equations

Recall from linear algebra that any system of linear equations a₁₁ x₁ + a₁₂ x₂ + a₁₃ x₃ = b₁
a₂₁ x₁ + a₂₂ x₂ + a₂₃ x₃ = b₂ can be written as a matrix equation or in even more reduced form as A x = b.

In this example, A is a 2 x 3 matrix, while x and b are column vectors (or 3 x 1 matrices). Number of rows first, then number of columns.

Theoretically from linear algebra we know the following. The first step is to understand the homogeneous equation, then we use that to understand the general situation.

• Homogeneous systems:
  • The system is homogeneous if b=0, i.e., b_k = 0 for all k.
  • A homogeneous system always has the "trivial" solution x=0.
  • If v is any other solution of a homogeneous system, then any multiple of v is a solution, and the sum of any two solutions is a solution.
  • Every matrix has a rank k, the number of linearly independent rows (or the number of linearly independent columns, which is always the same). Thus the rank of an m x n matrix is less than or equal to m and also to n (but may be strictly less).
    For example, the matrix has rank one, and the matrix has rank two, and no matrix of this form can have rank three.
  • The space of solutions of an m x n matrix has dimension n-rank(A). If rank(A) = n (i.e., the matrix has maximal rank), then the solution is unique.
  • The matrix is square if m=n. In this case the determinant of the matrix is a number det(A) (which we will worry about later) and the rank of the matrix is maximal if and only if det(A) ≠ 0.
• Inhomogeneous systems:
  
  
  • The system is inhomogeneous if b≠0, i.e., at least one b_k ≠ 0.
  • An inhomogeneous system may have no solution at all. In this case it is called inconsistent. For example, x+y=2, x+y=3 is inconsistent.
  • If an inhomogeneous system has one particular solution, then every other solution is the sum of this particular solution and a solution of the homogeneous equation.

This is theoretical. Notice what can go wrong numerically.

• The determinant may be very close to zero but not zero; numerically the computer may think it is zero. Or the opposite can happen.
• Alternatively the rows may be very close to dependent but not dependent. Numerically the computer may think the rows are dependent.

So the computer may find solutions that don't actually exist, or may be unable to find the solutions that do exist, due to roundoff errors. This is one major issue we will address in this chapter.

The other main issue is that even if such things don't happen, the number of computations required to solve the system may be very large. It's common to have a 100x100 matrix, for example, which may require a million computations to solve. So we want to make the algorithms as efficient as possible.

---

If you don't know anything else about the matrix, the simplest algorithm is Gaussian elimination (with backward substitution).

We subtract off successive multiples of previous rows to clear out all lower-triangular terms. We are left with a general upper triangular matrix. We do not reduce the matrix to the identity or reduced row echelon form.

Given a[i,j] and b[i], solve for x[i].

% First make the matrix upper-triangular.
for k from 1 to n-1 do
   for i from k+1 to n do
      m = a[i,k] / a[k,k]
      for j from k+1 to n do
         a[i,j] = a[i,j] - m*a[k,j]
      end do
      b[i] = b[i] - m*b[k]
   end do
end do

% Now solve for the x[i].
x[n] = b[n] / a[n,n]
for i from 1 to n-1 do
   x[n-i] = b[n-i]
   for j from n-i+1 to n do
      x[n-i] = x[n-i] - a[n-i,j]*x[j]
   end do
   x[n-i] = x[n-i] / a[n-i,n-i]
end do

Of course, in the second set of loops, one can do a backward "for" loop if one's program allows it (which all of our software does).

Let's count the number of multiplications and divisions when solving for the matrix:

• For each j-step, one multiplication. Total in the j-loop: (n-k).
• For each i-step, (n-k) multiplications from the j-loop, plus two more. So (n-k+2) for each i-step, for a total of (n-k+2)(n-k).
• Then since we are doing this from k=1 to k=n-1, we add up: Σ_k=1^n-1 (n-k)(n-k+2) = (2n³ + 3n² - 5n)/6.

Similarly we can count the number of additions and subtractions in solving for the matrix, which is Σ_k=1^n-1 (n-k)(n-k+1) = (n³ - n)/3.
We also want to count the number of multiplications and divisions in solving for x: Σ_i=1^n-1 (i + n) = (n²+n)/2.
And the number of additions and subtractions in solving for x: Σ_i=1^n-1 i = (n²-n)/2.
Observe:
The number of computations in solving for x is smaller than the number in solving for the matrix:
O(n²) to solve for x, vs. O(n³) to solve for the matrix.

This is easy to understand, since there are three loops in the matrix computation and only two loops in the x computation.

Roughly the number of operations is n³/3.

---

If we try to eliminate also the upper triangular terms, we get the Gauss-Jordan method:

1. Apply Gaussian elimination to reduce the matrix to upper triangular.
2. Once you have the upper triangular matrix, apply the same row reduction algorithm backwards to eliminate the terms above the diagonal. (Note there are fewer computations than in the first step since some rows are known to have zeroes in both places.
3. Stop when your matrix is diagonal.

You can write the same type of algorithm as above for this (and you will, for homework) and analyze the number of operations the same way.