MATH/APPM 4650
Class notes, April 23
Section 6.4 and 6.5: Determinants and factorization

Determinants

The determinant of a 2x2 matrix is

| a11	a12 |
| a21	a22 |

= a₁₁a₂₂ - a₁₂a₂₁

The determinant of a 3x3 matrix is

| a11	a12	a13 |
| a21	a22	a23 |
| a31	a32	a33 |

= a₁₁a₂₂a₃₃ + a₁₂a₂₃a₃₁ + a₁₃a₂₁a₃₂ - a₁₂a₂₁a₃₃ - a₁₃a₂₂a₃₁ - a₁₁a₂₃a₃₂

In general, the determinant of an nxn matrix has n! terms (the number of permutations of n elements) and is given by Σ_Sn sgn(σ) a_1σ(1)a_2σ(2)...a_nσ(n) where σ is a permutation of {1,2,...,n} and sgn(σ) is its sign (positive if the number of transpositions is even, negative if it's odd).

Obviously to compute the determinant from the definition takes roughly n! computations. (n or n-1 multiplications for each summand, then n!-1 additions.)

10! = 3,628,800, so computing the determinant is extremely time-consuming.

Thus for example, one would almost never want to compute the determinant before trying to invert a matrix. It's much faster to just try to invert it first and stop if that process fails.

Of course, there are faster ways to compute the determinant.

Expanding on a row (or a column) is not faster.
To compute the determinant of a matrix of size n, you need to do n multiplications and n determinants of size n-1.
Each of those n-1 determinants requires n-1 multiplications and the computation of n-1 determinants of size n-2...
And so on.
We still get an algorithm of size n! or so.

Instead, we can use the Gaussian elimination algorithm (which we know takes O(n³) operations) to reduce the matrix to an upper-triangular matrix.

The determinant does not change when any of the Gaussian elimination operations are used (except switching rows, which multiplies the determinant by -1 each time).

The determinant of an upper triangular matrix is the product of the diagonal terms. So after a reduction, it's only n multiplications to get the determinant.

This is by far the most practical way to compute the determinant of a matrix (unless you have some other information about the matrix).

---

Factorization

Recall that it is easy (O(n²) computations) to solve Ux = b for x when U is an upper-triangular matrix.

% Solve Ux=b for x, when U is an upper triangular matrix with no zeros on the diagonal.
x[n] = b[n] / U[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] - U[n-i,j]*x[j]
   end do
   x[n-i] = x[n-i] / U[n-i,n-i]
end do

Similarly it is easy to solve Lx = b for x if L is a lower-triangular matrix.

% Solve Lx=b for x, when L is a lower triangular matrix with no zeros on the diagonal.
x[1] = b[1] / L[1,1]
for i from 1 to n-1 do
   x[i] = b[i]
   for j from 1 to i-1 do
      x[i] = x[i] - L[i,j]*x[j]
   end do
   x[i] = x[i] / L[i,i]
end do

In general, most matrices can be written in the form A = LU, for some lower-triangular matrix L and upper-triangular matrix U.

So if we want to solve Ax = b, we can write LUx = b.

Letting y = Ux, we can think of this as two matrix equations: Ly = b and Ux = y.

First solve Ly = b for y. Then solve Ux = y for x.

Both solutions take O(n²) operations, so we have solved a matrix using only O(n²) operations.

Of course, the problem is how to find L and U given A. This takes O(n³) operations.

No free lunch!



If we look at exactly how we did the row reduction process, we can think of each step of Gaussian elimination as being a matrix multiplication.

1. "Multiply the second row of A by 5" is the same as "Replace A with M₁A, where M₁ is the following matrix."
   

   / 1	0 \
   \ 0	5 /

2. "Add 7 times row one to row two" is the same as "Replace A with M₂A, where M₂ is the following matrix."
   

   / 1	0 \
   \ 7	1 /

So for example, with a 3x3 matrix, to obtain the upper triangular form, we would do three of the second kind of operation.
U = M₃M₂M₁A where

	/   1	0	0 \
M1   =	| -m12	1	0 |
	\   0	0	1 /

	/     1	0	0 \
M2   =	|      0	1	0 |
	\ -m13	0	1 /

	/ 1	0	0 \
M3   =	| 0	1	0 |
	\ 0	-m23	1 /

(Important: notice that m₂₃ is the multiplier for the reduced matrix; that is, you clear out the first column first, then compute a₃₂/a₂₂ using the new values of a.

This means that A = M₁^-1M₂^-1M₃^-1U

But notice that

	/ 1	0	0 \
M1-1M2-1M3-1    =	| m12	1	0 |
	\ m13	m23	1 /

So we have factored A = LU, where L is a lower-triangular matrix (with 1s on the diagonal) and U is an upper-triangular matrix.

This can always be done, as long as we never have to switch rows.

Switching rows is a matrix like

/ 0	1 \
\ 1	0 /

and this messes up the lower-triangular property.

So the algorithm is:
% Factor a matrix into LU form.
L[1,1] = 1
for k from 2 to n do
   L[k,k] = 1
   for i from k to n do
      L[i,k-1] = a[i,k-1] / a[k-1,k-1]
      for j from k to n do
         a[i,j] = a[i,j] - L[i,k-1]*a[k-1,j]
      end do
   end do
end do
% Save the new entries of A as U.

Why are we doing this?

If we want to solve the equation Ax=b multiple times, for the same A but different b, we should not use the full Gaussian elimination algorithm each time.
Do the LU factorization once, then solve Ly=b, Ux=y repeatedly.