Mathematics 372 – Numerical Linear Algebra Midterm Problem Set Solutions March 23, 2007

I. (More on Gaussian elimination and LU factorizations).

A) (10) Show that an n × n matrix A has an LU factorization (equivalently, Gaussian elimination is possible with no row interchanges) if and only if for each 1 ≤ k ≤ n, the upper left corner k × k submatrix

A_{k} = ( \begin{matrix}a_{11} & ⋯ & a_{1k} \\ ⋮ & & ⋮ \\ a_{k1} & ⋯ & a_{kk} \end{matrix} )

is nonsingular.

Solution: ⇒: First assume that Gaussian elimination on A is possible with no row interchanges. We will argue by induction on k. The base case is k = 1. If we do not have to interchange rows to find a pivot element in column 1, then a_{11} = det(A_{1}) ≠ 0. So the base case is established. Now assume that we have shown that det(A_{ℓ}) ≠ 0 for ℓ < k. We do the row operations (without row interchanges) to create the zeroes below the main diagonals in columns 1 through ℓ. Since each of these row operations is of type 1 (replacing a row by that row plus a multiple of another row), the determinants of the matrix A and of the principal minors are not changed. In particular, the determinant of the minor A_{ℓ+1} will be the same as the determinant of the “eliminated” (ℓ + 1) × (ℓ + 1) matrix:

( \begin{matrix}a_{11} & a_{12} & a′_{13} & ⋯ & a′_{1ℓ} & a_{1,ℓ+1} \\ 0 & a′_{22} & a′_{23} & ⋯ & a′_{2ℓ} & a′_{2,ℓ+1} \\ 0 & 0 & a′_{33} & ⋯ & a′_{3ℓ} & a′_{3,ℓ+1} \\ ⋮ & ⋮ & ⋮ & ... & ⋮ & ⋮ \\ 0 & 0 & 0 & ⋯ & a′_{ℓℓ} & a′_{ℓ,ℓ+1} \\ 0 & 0 & 0 & ⋯ & 0 & a′_{ℓ+1,ℓ+1} \end{matrix} )

The induction hypothesis is that

det(A_{ℓ}) = a_{11}a′_{22} ⋯ a′_{ℓℓ} ≠ 0.

If we can use the diagonal entry as pivot in column ℓ + 1, then a_{ℓ+1,ℓ+1} ≠ 0. Hence

det(A_{ℓ+1}) = det(A_{ℓ})a′_{ℓ+1,ℓ+1} ≠ 0

also.

⇐: For the converse, we essentially just argue in the reverse direction. If all the det(A_{k}) ≠ 0 for 1 ≤ k ≤ n, then doing the reductions as above, after completing each column, the diagonal entry in the next column must be nonzero. Hence we never have to interchange rows in the Gaussian elimination.

B) (10) Suppose that A is a matrix that satisfies the condition in part A, and for any k, 1 ≤ k ≤ n, partition A as

A = ( \begin{matrix}A_{k} & B \\ C & D \end{matrix} ) ,

where B is k × (n − k), C is (n − k) × k, and D is (n − k) × (n − k). Show that there is a unique (n − k) × k matrix M such that

( \begin{matrix}I_{k} & 0 \\ −M & I_{n−k} \end{matrix} ) ( \begin{matrix}A_{k} & B \\ C & D \end{matrix} ) = ( \begin{matrix}A_{k} & B̂ \\ 0 & D̂ \end{matrix} )

for some matrices B̂ and D̂. Also show that D̂ = D − CA_{k}^{−1}B.

Solution: Matrices with these block decompositions can be multiplied using the same formulas as for 2 × 2 scalar matrices. In particular the entry in the “row 2, column 1” of the product on the left side of the equation to be proved is

−MA_{k} + I_{n−k}C = −MA_{k} + C.

Because A_{k} is assumed to be invertible, there is only one matrix that makes this equation true, namely:

(∗)  M = CA_{k}^{−1}.

The entry in “row 1, column 2” of the product is B̂ = B. From (*), −MB + I_{n−k}D = D̂, so D̂ = D − MB = D − CA_{k}^{−1}B as desired.

C) (5) Deduce that A has the block LU decomposition

A = ( \begin{matrix}I_{k} & 0 \\ M & I_{n−k} \end{matrix} ) ( \begin{matrix}A_{k} & B̂ \\ 0 & D̂ \end{matrix} ) .

Solution: This follows because

( \begin{matrix}I_{k} & 0 \\ −M & I_{n−k} \end{matrix} )^{−1} = ( \begin{matrix}I_{k} & 0 \\ M & I_{n−k} \end{matrix} )

(check by direct matrix multiplication).

D) (5) By definition, the matrix A_{k} satisfies the condition in part A. It can be shown that D̂ does as well. Hence each of these matrices also has an LU-decomposition:

A_{k} = L_{1}U_{1}, and D̂ = L_{2}U_{2}. Express the LU-decomposition of A in terms of L_{1}, L_{2}, M, U_{1}, U_{2}, and B.

Solution: We have, by regrouping terms using associativity of matrix multiplication,

A = ( \begin{matrix}I_{k} & 0 \\ M & I_{n−k} \end{matrix} ) ( \begin{matrix}A_{k} & B \\ 0 & D̂ \end{matrix} )

= ( \begin{matrix}I_{k} & 0 \\ M & I_{n−k} \end{matrix} ) ( \begin{matrix}L_{1}U_{1} & B \\ 0 & L_{2}U_{2} \end{matrix} )

= ( \begin{matrix}I_{k} & 0 \\ M & I_{n−k} \end{matrix} ) ( \begin{matrix}L_{1} & 0 \\ 0 & L_{2} \end{matrix} ) ( \begin{matrix}U_{1} & L_{1}^{−1}B \\ 0 & U_{2} \end{matrix} )

= ( \begin{matrix}L_{1} & 0 \\ ML_{1} & L_{2} \end{matrix} ) ( \begin{matrix}U_{1} & L_{1}^{−1}B \\ 0 & U_{2} \end{matrix} ) .

This is an ordinary LU factorization for A.

II. (The matrix 2-norm, revisited.) By the general theory of the matrix operator norms associated to vector norms on R^{n}, the matrix 2-norm of A ∈ M_{n×n}(R) is given by

(1)  ∥A∥_{2} = max_{∥u∥_{2}=1} ∥Au∥_{2}.

Your goal in this problem is to prove the formula

(2)  ∥A∥_{2} = \sqrt{max\{|λ| : λ is an eigenvalue of A^{t}A\}}

that we mentioned in class. The equation (1) shows that we need to solve a constrained optimization problem to determine ∥A∥_{2}.

A) (5) What standard method from multivariable calculus applies to solve constrained optimization problems of the following form:

maximize : f(x_{1}, …, x_{n})
subject to : g(x_{1}, …, x_{n}) = 0?

Describe how the method works in general.

Solution: The method is the Lagrange multiplier method for constrained optimization. One way to say the process is that you consider the system of n + 1 equations in n + 1 variables x_{1}, …, x_{n}, λ obtained from the components of the vector equation ∇f = λ∇g, together with the gradient equation. The maximum on the constraint set will be among the solutions of this system.

B) (15) For our purposes in (1), it will be most convenient to square both the function ∥Au∥_{2} we seek to maximize and the constraint equation ∥u∥_{2} = 1. Writing u = (u_{1}, …, u_{n})^{t}, apply the method with f(u_{1}, …, u_{n}) = ∥Au∥_{2}^{2} and g(u_{1}, …, u_{n}) = u_{1}^{2} + ⋯ + u_{n}^{2} − 1. Show that the equations for the constrained maximum imply that the vector u where f achieves its maximum must be an eigenvector of A^{t}A, and that the multiplier λ in the method must be the corresponding eigenvalue.

Solution: With

f = ∥Au∥_{2}^{2} = ⟨Au, Au⟩ = ⟨u, A^{t}Au⟩

and g = ∥u∥_{2}^{2} = ⟨u, u⟩, in the Lagrange equations from ∇f = λ∇g we have the equations

\frac{∂f}{∂u_{i}} = λ\frac{∂g}{∂u_{i}}

for i = 1, …, n. From the form of g, the right side is easy: 2λu_{i}. The left side is somewhat messy, but simplest if you use the form f = ⟨u, A^{t}Au⟩. When we take this expression and compute the derivative with respect to u_{i}, we get 2 times the ith component of A^{t}Au. Hence the equation ∇f = λ∇g is just 2A^{t}Au = 2λu, or A^{t}Au = λu. Since u is a unit vector from the constraint equation, it is nonzero. Hence u is an eigenvector of A^{t}A with eigenvalue λ.

C) (5) Deduce that (2) holds, giving the formula for ∥A∥_{2}.

Solution: From the previous part, the solutions of the Lagrange equations will consist of (u, λ) where u is a unit eigenvector of A^{t}A, and λ is the corresponding eigenvector. The maximum is obtained for the one of these that gives the largest value for f. But notice,

f(u) = ⟨Au, Au⟩ = ⟨u, A^{t}Au⟩ = ⟨u, λu⟩ = λ⟨u, u⟩ = λ.

In other words, the value of f at the eigenvector u is just the corresponding eigenvalue of A^{t}A. Hence the maximum value of f will be

max{λ : λ is an eigenvalue of A^{t}A}

(note the absolute value is actually unnecessary because A^{t}A is symmetric and positive semidefinite, so all the eigenvalues are reals with λ ≥ 0). The desired equation (2) comes from taking square roots:

∥A∥_{2} = \sqrt{max\{λ : λ is an eigenvalue of A^{t}A\}}.

D) (10) Recall the 60 × 60 adjacency matrix B of the “bucky ball” from Lab/Problem Set 3. Use (2) above and MATLAB to compute ∥B∥_{2}. Some technical notes: (a) The MATLAB command eig returns a list of all the eigenvalues of a square matrix; (b) the built-in matrix bucky is a sparse matrix. If you want to check your answer using the norm command, you will need to convert bucky to “full” form.

Solution: See MATLAB solutions.

E) (10) If A = (a_{ij}), the matrix Frobenius norm is defined by

∥A∥_{F} = \sqrt{∑_{i=1}^{n} ∑_{j=1}^{n} |a_{ij}|^{2}}

(note: the absolute values here are unnecessary for real matrices, but necessary if some of the entries in A are complex numbers with nonzero imaginary parts). Show that for all matrices A

(3)  ∥A∥_{2} ≤ ∥A∥_{F} ≤ \sqrt{n}∥A∥_{2}

Solution: For the first inequality, subdivide A into its rows and call these vectors A_{1}, …, A_{n}. Then for some unit vector u,

∥A∥_{2}^{2} = ∥Au∥_{2}^{2}
= ⟨A_{1}, u⟩^{2} + ⋯ + ⟨A_{n}, u⟩^{2}
≤ ∥A_{1}∥_{2}^{2}∥u∥_{2}^{2} + ⋯ + ∥A_{n}∥_{2}^{2}∥u∥_{2}^{2} (Cauchy − Schwarz)
= ∥A_{1}∥_{2}^{2} + ⋯ + ∥A_{n}∥_{2}^{2} (since ∥u∥_{2} = 1)

This last line is the sum of the squares of all the entries in the matrix, which is ∥A∥_{F}^{2} . Hence ∥A∥_{2} ≤ ∥A∥_{F} for all n × n matrices A. For the second inequality, split A into its columns, A^{(1)}, …, A^{(n)}. Then the Frobenius norm squared is the sum of the squares of the 2-norms of the columns. Using the fact that A^{(j)} = Ae_{j} and standard properties of matrix norms,

∥A∥_{F}^{2} = ∥A^{(1)}∥_{2}^{2} + … + ∥A^{(n)}∥_{2}^{2}
= ∥Ae_{1}∥_{2}^{2} + … + ∥Ae_{n}∥_{2}^{2}
≤ ∥A∥_{2}∥e_{1}∥_{2}^{2} + … + ∥A∥_{2}∥e_{1}∥_{2}^{2}
≤ n∥A∥_{2}

(since ∥e_{i}∥_{2}^{2} = 1 for all of the standard basis vectors e_{i}). Hence taking square roots, ∥A∥_{F} ≤ \sqrt{n}∥A∥_{2}.

F) (5) Compute the Frobenius norm of the bucky ball matrix B using MATLAB and verify that (3) holds in this case.

Solution: See MATLAB solutions.

III. (Roundoff considerations and least squares.) In many situations, the way a problem is set up can have a significant effect on the susceptibility to roundoff errors when we apply standard methods. But fortunately, there are often ways to minimize those effects. For instance, suppose that we want to compute the best-fit line for the following data points

\begin{matrix}x_{i} & 1.001 & 1.002 & 1.003 & 1.004 & 1.005 & 1.006 & 1.007 \\ y_{i} & 3.664 & 3.789 & 3.891 & 4.022 & 4.233 & 5.200 & 5.329 \end{matrix}

by solving the normal equations X^{t}X ( \begin{matrix}m \\ b \end{matrix} ) = X^{t}Y for the least squares problem.

A) (5) Using MATLAB, find the ∞-norm condition number of the coefficient matrix X^{t}X in the system of normal equations. What does your answer indicate about this system of equations?

Solution: See MATLAB solutions.

B) (10) Our standard way of setting up the normal equations corresponds to using the basis {x, 1} for the vector space of linear polynomials. But we could use any other basis we like too. Suppose we use {300(x − 1.004), 1} instead (that is, translate and rescale the x_{i} values). What is the ∞-norm condition number of the coefficient matrix in the normal equations now? Show that the best-fit line can still be computed if we set the problem up this way.

Solution: See MATLAB solutions.

C) (5) Explain your results in A and B geometrically.

Solution: The idea is that the matrix X^{t}X is nearly singular (this is what the large condition numbers mean. The lines in R^{2} coming from the two equations in the system

X^{t}X ( \begin{matrix}m \\ b \end{matrix} ) = X^{t}Y are nearly parallel. By rescaling (multiplying the x-coordinates

by 300), the slopes are made much smaller, so perturbing the lines shifts the position of the intersection less in a relative sense. This is indicated by the smaller condition number of the second matrix.