Matrix Algorithms (15)

Michael Littman November 17th, 1998

1 MATRIX APPLICATIONS

1.1 Scientific Computing

Use computers to simulate the real world:

• Semiconductor process

• Aircraft/automobile design

• Electrocardiography (heart beats) How?

1. Set up the mathematical model.

2. Parse the mathematical model into a system of linear equations.

3. Solve the system of linear equations.

1.2 Learning from Data

A common task in statistics, data mining, machine learning, etc. is to find a mathematical function that fits a set of data. For example: Given temperature, wind speed, humidity, day of the year, etc., predict the strength of a hurricane. Or: Given n, predict the running time of selection sort:

This shows the “best” linear function (0.001779818n − 3.596).

2 MATRICES AND VECTORS

2.1 Definitions

Matrix:

• A 2-dimensional array of real numbers

• An m by n matrix has m rows and n columns

A_{m×n} = [ \begin{matrix}a_{11} & a_{12} & ⋯ & a_{1n} \\ a_{21} & a_{22} & ⋯ & a_{2n} \\ ⋮ & ⋮ & & ⋮ \\ a_{m1} & a_{m2} & ⋯ & a_{mn} \end{matrix} ]

• A[i, j] = a_{ij}, element at ith row and jth column. Vector:

• A 1-dimensional array of real numbers

• A row vector (of size n) — a matrix with 1 row

u^{T} = [u_{1} u_{2} ⋯ u_{n}]

• A column vector (of size m) — a matrix with 1 column

v = [ v_{1} ; v_{2} ; ⋮ ; v_{m} ]

• By default, a vector is a column vector

2.2 Some Special Matrices

Square matrix: Same number of rows and columns (n × n). Identity matrix, I_{n×n}:

I[i, j] = \{ \begin{matrix}1 & if & i = j \\ 0 & otherwise \end{matrix}

Diagonal matrix, D_{n×n}

D[i, j] = 0 if i ≠ j

Tridiagonal matrices, T_{n×n}

T [i, j] = 0 if |i − j| > 1

Lower-triangular matrices, L_{n×n}

L[i, j] = 0 if i < j

Upper-triangular matrices, U_{n×n}

U[i, j] = 0 if i > j

Pictures...

2.3 Basic Matrix Operations

Transpose: A^{T} • Algorithm: A^{T} [i, j] = A[j, i], exchange the row and column indices.

• # operations = Θ(mn) Scalar multiplication: c · A_{m×n}

• Precondition: c is a real number.

• Algorithm: (c · A)[i, j] = c · A[i, j], multiply each element of A by c.

• # operations = Θ(mn) Addition: A_{m×n} + B_{m×n}

• Precondition: A and B have the same shape.

• Algorithm: (A + B)[i, j] = A[i, j] + B[i, j], add element by element.

• # operations = Θ(mn)

Dot product (of vectors): u^{T} v

• Preconditions:

– A row vector (u^{T}) times a column vector (v);

– Two vectors have the same size (n).

• Algorithm: u^{T} v = u_{1}v_{1} + u_{2}v_{2} + ⋯ + u_{n}v_{n} ≡ ∑_{i=1}^{n} u_{i}v_{i}

• # operations = Θ(n) Norm (l^{2}-norm) of a vector: ∥u∥ = \sqrt{u^{T} u}. How can we express matrix subtraction in terms of these operations?

A − B = A + ((−1) · B).

2.4 Matrix Multiplication

A_{m×n} · B_{n×p}:

• Precondition: # columns of A equals # rows of B.

• Algorithm: (A · B)[i, j] = A[i, 1]B[1, j] + A[i, 2]B[2, j] + ⋯ + A[i, n]B[n, j]

• Alternatively, write A as an array of row vectors

A = [ a_{1}^{T} ; a_{2}^{T} ; ⋮ ; a_{m}^{T} ]

and B as an array of column vectors

B = [b_{1} b_{2} ⋯ b_{p}] ,

then

A · B = [ \begin{matrix}a_{1}^{T} b_{1} & a_{1}^{T} b_{2} & ⋯ & a_{1}^{T} b_{p} \\ a_{2}^{T} b_{1} & a_{2}^{T} b_{2} & ⋯ & a_{2}^{T} b_{p} \\ ⋮ & ⋮ & & ⋮ \\ a_{m}^{T}b_{1} & a_{m}^{T}b_{2} & ⋯ & a_{m}^{T}b_{p} \end{matrix} ] .

• # operations = Θ(mnp):

– A · B involves (m × p) dot products, and – Each dot product needs O(n) operations.

2.5 Matrix Algebra

For all matrices A_{m×n}, B_{m×n} and scalars c:

• c(A + B) = cA + cB

• A + B = B + A

• (A^{T})^{T} = A

• A − A = 0

• (A + B)^{T} = A^{T} + B^{T} In addition, for all matrices C_{p×m}:

• (AB)^{T} = B^{T} A^{T}

• C(A + B) = CA + CB However, AB ≠ BA in general (even if both are square):

[ \begin{matrix}1 & 2 \\ 3 & 4 \end{matrix} ] [ \begin{matrix}4 & 2 \\ 3 & 1 \end{matrix} ] = [ \begin{matrix}10 & 4 \\ 24 & 10 \end{matrix} ] ,

but

[ \begin{matrix}4 & 2 \\ 3 & 1 \end{matrix} ] [ \begin{matrix}1 & 2 \\ 3 & 4 \end{matrix} ] = [ \begin{matrix}10 & 16 \\ 6 & 10 \end{matrix} ] .

2.6 Solving Equations

Say A and C are m × n and b and x are n-vectors. Solve for vector x:

(C^{T} A)^{T} b + x = b − 2x
A^{T} Cb + x = b − 2x
x + 2x = b − A^{T} Cb
3x = b − A^{T} Cb
x = 1/3(b − A^{T} Cb)

Can solve lots of equations this way. What about Ax = b? Need a notion of matrix inverse.

3 SOLVING LINEAR SYSTEMS

3.1 Gaussian Elimination

Basic form: A · x = b. An example:

A = [ \begin{matrix}6 & −2 & 2 & 4 \\ 12 & −8 & 6 & 10 \\ 3 & −13 & 9 & 3 \\ −6 & 4 & 1 & −18 \end{matrix} ] ,

x = [ \begin{matrix}x_{1} \\ x_{2} \\ x_{3} \\ x_{4} \end{matrix} ] ,

b = [ \begin{matrix}12 \\ 34 \\ 27 \\ −38 \end{matrix} ] .

\{ \begin{matrix}6x_{1} & − & 2x_{2} & + & 2x_{3} & + & 4x_{4} \\ 12x_{1} & − & 8x_{2} & + & 6x_{3} & + & 10x_{4} \\ 3x_{1} & − & 13x_{2} & + & 9x_{3} & + & 3x_{4} \\ − 6x_{1} & + & 4x_{2} & + & 1x_{3} & − & 18x_{4} \end{matrix}  = = = = \begin{matrix}12 \\ 34 \\ 27 \\ −38 \end{matrix}

1. Eliminate x_{1} in equations 2,3 and 4

\{ \begin{matrix}6x_{1} & − & 2x_{2} & + & 2x_{3} & + & 4x_{4} \\ & − & 4x_{2} & + & 2x_{3} & + & 2x_{4} \\ & − & 12x_{2} & + & 8x_{3} & + & x_{4} \\ & & 2x_{2} & + & 3x_{3} & − & 14x_{4} \end{matrix}  = = = = \begin{matrix}12 \\ 10 \\ 21 \\ −26 \end{matrix}

2. Eliminate x_{2} in equations 3 and 4

\{ \begin{matrix}6x_{1} & − & 2x_{2} & + & 2x_{3} & + & 4x_{4} \\ & − & 4x_{2} & + & 2x_{3} & + & 2x_{4} \\ & & & & 2x_{3} & − & 5x_{4} \\ & & & & 4x_{3} & − & 13x_{4} \end{matrix}  = = = = \begin{matrix}12 \\ 10 \\ −9 \\ −21 \end{matrix}

3. Eliminate x_{3} in equation 4

\{ \begin{matrix}6x_{1} & − & 2x_{2} & + & 2x_{3} & + & 4x_{4} \\ & − & 4x_{2} & + & 2x_{3} & + & 2x_{4} \\ & & & & 2x_{3} & − & 5x_{4} \\ & & & & & − & 3x_{4} \end{matrix}  = = = = \begin{matrix}12 \\ 10 \\ −9 \\ −3 \end{matrix}

4. Solve the above system starting from the last equation

\{ \begin{matrix}x_{1} \\ x_{2} \\ x_{3} \\ x_{4} \end{matrix}  = = = = \begin{matrix}1 \\ −3 \\ −2 \\ 1 \end{matrix}

3.2 Expressing the Algorithm

There is a more systematic way of describing the Gaussian elimination process.

• Easy-to-solve systems – A is a lower-triangular matrix

[ \begin{matrix}a_{11} & 0 & 0 & ⋯ & 0 \\ a_{21} & a_{22} & 0 & ⋯ & 0 \\ a_{31} & a_{32} & a_{33} & ⋯ & 0 \\ ⋮ & ⋮ & ⋮ & & ⋮ \\ a_{n1} & a_{n2} & a_{n3} & ⋯ & a_{nn} \end{matrix} ] · [ x_{1} ; x_{2} ; x_{3} ; ⋮ ; x_{n} ] = [ b_{1} ; b_{2} ; b_{3} ; ⋮ ; b_{n} ] .

It can be solved as

x_{1} = b_{1}/a_{11},
x_{2} = (b_{2} − a_{21}x_{1})/a_{22},
x_{3} = (b_{3} − a_{31}x_{1} − a_{32}x_{2})/a_{33},
⋮ ⋮
x_{n} = (b_{n} − a_{n1}x_{1} − a_{n2}x_{2} − ⋯ − a_{n n−1}x_{n−1})/a_{nn}.

This is called the forward substitution.

– A is an upper-triangular matrix

[ \begin{matrix}a_{11} & a_{12} & a_{13} & ⋯ & a_{1n} \\ 0 & a_{22} & a_{23} & ⋯ & a_{2n} \\ 0 & 0 & a_{33} & ⋯ & a_{3n} \\ ⋮ & ⋮ & ⋮ & & ⋮ \\ 0 & 0 & 0 & ⋯ & a_{nn} \end{matrix} ] · [ x_{1} ; x_{2} ; x_{3} ; ⋮ ; x_{n} ] = [ b_{1} ; b_{2} ; b_{3} ; ⋮ ; b_{n} ] .

It can be solved as

x_{n} = b_{n}/a_{nn},
x_{n−1} = (b_{n−1} − a_{n−1,n}x_{n})/a_{n−1,n−1},
x_{n−2} = (b_{n−2} − a_{n−2,n}x_{n} − a_{n−2,n−1}x_{n−1})/a_{n−2,n−2},
⋮ ⋮
x_{1} = (b_{1} − a_{1n}x_{n} − a_{1,n−1}x_{n−1} − a_{1,n−2}x_{n−2} ⋯ − a_{12}x_{2})/a_{11}.

This is called the back substitution – Constraint: a_{11}, a_{22}, …, a_{nn} ≠ 0

– # operations for both cases — Θ(n^{2}).

• Strategy If we can write A = LU, with L lower-triangular and U upper-triangular, then what we need to do is

1. Solve z in Lz = b, then

2. Solve x in Ux = z

3.3 LU-factorization

• Algorithm 1. If a_{11} = 0, fail!

2. If n = 1, return L = 1 and U = a_{11}.

3. Write A as

A = [ \begin{matrix}a_{11} & a_{12} & ⋯ & a_{1n} \\ a_{21} & a_{22} & ⋯ & a_{2n} \\ ⋮ & ⋮ & & ⋮ \\ a_{n1} & a_{n2} & ⋯ & a_{nn} \end{matrix} ] = [ \begin{matrix}a_{11} & r^{T} \\ c & A′ \end{matrix} ] ,

4. Define A_{1} = A′ − \frac{c · r^{T}}{a_{11}} . (Schur complement, see CLR p.755)

5. Recursively compute the LU-factorization of A_{1}

A_{1} = L_{1}U_{1}

6. Return

L = [ \begin{matrix}1 & 0 \\ c/a_{11} & L_{1} \end{matrix} ] , U = [ \begin{matrix}a_{11} & r^{T} \\ 0 & U_{1} \end{matrix} ]

• Correctness of the algorithm

– Show L · U = A 1. Base case — L · U = 1 · a_{11} = a_{11} = A

2. Otherwise —

L · U = [ \begin{matrix}1 & 0 \\ c/a_{11} & L_{1} \end{matrix} ] · [ \begin{matrix}a_{11} & r^{T} \\ 0 & U_{1} \end{matrix} ]
= [ \begin{matrix}a_{11} & r^{T} \\ c & L_{1} · U_{1} + c · r^{T} /a_{11} \end{matrix} ]

= [ \begin{matrix}a_{11} & r^{T} \\ c & A_{1} + c · r^{T} /a_{11} \end{matrix} ]
= [ \begin{matrix}a_{11} & r^{T} \\ c & A′ \end{matrix} ]
= A

– Show L is lower-triangular and U is upper-triangular (Exercise)

– Termination: # recursive calls = n−1

• Operation count – Computing A_{1} (Schur complement) costs Θ(n^{2}) operations

– T (n) = T (n − 1) + Θ(n^{2})

– Solution?

T (n) = Θ(n^{3}).

3.4 P LU-factorization

P LU is LU with pivoting (although “P” stands for permutation).

• Permutation matrices P_{n×n}:

– Each row or each column has exactly one 1 and 0 for the rest.

– Example: [ \begin{matrix}0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{matrix} ]

– Property: P^{T} P = P P^{T} = I

• What if a_{11} = 0 happens in LU-factorization? – Want to find a permutation matrix P, such that A = P LU.

• Algorithm 1. If n = 1, return P = 1, L = 1, U = a_{11}.

2. Find a non-zero element a_{i1} in the first row of A. (If not, then the system is not solvable.)

3. Define a permutation matrix Q_{n×n} such that – q_{1i} = q_{i1} = 1,

– q_{jj} = 1 for j ≠ 1, j ≠ i,

– All the other elements of Q are 0.

4. Let B = QA. Notice that b_{11} = a_{i1} ≠ 0.

5. Write B as

B = [ \begin{matrix}b_{11} & r^{T} \\ c & B′ \end{matrix} ] .

6. Define B_{1} = B′ − c · r^{T} /b_{11}. Recursively compute B_{1} = P_{1}L_{1}U_{1} 7. Return

P = Q^{T} [ \begin{matrix}1 & 0 \\ 0 & P_{1} \end{matrix} ] , L = [ \begin{matrix}1 & 0 \\ P_{1}^{T} c/b_{11} & L_{1} \end{matrix} ] , U = [ \begin{matrix}b_{11} & r^{T} \\ 0 & U_{1} \end{matrix} ]

• Correctness of the algorithm? By induction.

• Operation count? Same as LU.

3.5 Summary

For an n × n matrix A, the following are equivalent.

1. A has a P LU-factorization.

2. For any n-vector b, the system Ax = b has a solution.

3. A has an inverse.

4. A has non-zero determinant.

5. A has full rank. Read CLR Ch 31.1 for definitions of inverse, determinant and rank.

4 LEAST SQUARES

4.1 Overconstrained Systems

We discussed how to solve Ax = b, and it only works when A is square (same number of equations as variables). Consider the case of finding the best line to fit a set of points. We want to find x_{1} and x_{2} such that x_{1}n_{i} + x_{2} = t_{i} where n_{i} is problem size and t_{i} is time. We can’t solve this exactly (in general) if we have more than two data points (which is the only useful case anyway!).

4.2 Least Squares

So, what solution do we pick? A classic example is to find the x_{1} and x_{2} that give the best values in a “least squares” sense:

min_{x_{1},x_{2}} ∑_{i}(x_{1}n_{i} + x_{2} − t_{i})^{2}.

This scoring scheme is reasonable because it is optimized for a perfect fit, and the worse the predictions, the worse the score. Note: You can think of this as being a generalization of the notion of an average. Recall that min_{x_{1}} ∑_{i}(x_{1} − t_{i})^{2} = ∑_{i} t_{i}/n.

4.3 Least Squares in Matrices

We can also express the problem this way: Let A be a matrix with the n_{i} column as its first column and a column of all 1s for its second column. Let x be the column vector containing x_{1} and x_{2} and b be a column vector of the y_{i} values. We want to find the x that minimizes (Ax − b)^{2}. In this general form, we can also consider fitting other functions. For example, x_{1}n_{i}^{2} +x_{2}n_{i} + x_{3} ≈ t_{i} can be expressed by putting n_{i}^{2} values in the first column of A, n_{i} in the second column of A, and 1s in the third column.

4.4 Minimizing The Score

Recall from calculus that a differentiable function is minimized at the point at which the derivative is zero. Thus, we want to find

D((Ax − b)^{2}) = 0
2D(Ax − b)^{T} (Ax − b) = 0
A^{T} (Ax − b) = 0
A^{T} Ax = A^{T} b.

Cool! It’s just another system of linear equations to solve. (Thank goodness derivatives are linear operators!)

4.5 Example

A = [ \begin{matrix}1000 & 1 \\ 2000 & 1 \\ 3000 & 1 \\ 4000 & 1 \\ 5000 & 1 \\ 6000 & 1 \\ 7000 & 1 \\ 8000 & 1 \\ 9000 & 1 \\ 10000 & 1 \end{matrix} ] ,

b = [ \begin{matrix}0.16 \\ 0.63 \\ 1.44 \\ 2.56 \\ 3.99 \\ 5.73 \\ 7.84 \\ 10.28 \\ 12.98 \\ 16.32 \end{matrix} ] .

Form A^{T} Ax = A^{T} b, solve for x:

x = [ \begin{matrix}0.001779818 \\ −3.596 \end{matrix} ] .

Thus, t ≈ 0.001779818n − 3.596 is the best fit to the data in a least squares sense. The best quadratic is much better: t ≈ 1.67803 × 10^{−7}n^{2} − 6.601515 × 10^{−5}n + 0.09566667.