Mathematics 372 – Numerical Linear Algebra Solutions for Final Problem Set May 9, 2007

I. Some additional properties of matrix norms. Recall that from Theorem 4.2.1 in Watkins, using the SVD, we know ∥A∥_{2} = σ_{1}, the largest singular value of A. This gives additional properties and estimates for the matrix 2-norm.

A) (10) On the midterm problem set, recall that we showed

∥A∥_{F} ≤ \sqrt{n} · ∥A∥_{2}

for all n × n matrices. Show the following more general and sharper form of this inequality: For all A ∈ M_{n×m}(R),

∥A∥_{F} ≤ \sqrt{rank(A)} · ∥A∥_{2}.

Solution: Say rank(A) = r. Let A = UΣV^{t} be an svd for A. Then by the result of Exercise 4.2.3,

∥A∥_{F}^{2} = ∥Σ∥_{F}^{2} = σ_{1}^{2} + ⋯ + σ_{r}^{2} ≤ rσ_{1}^{2} = rank(A)σ_{1}^{2} = rank(A)∥A∥_{2}^{2}.

Taking square roots completes the proof.

B) (10) Show that for all matrices A ∈ M_{n×m}(R),

∥A∥_{2} ≤ \sqrt{∥A∥_{1}∥A∥_{∞}}

(this is sometimes useful for estimating ∥A∥_{2} without computing it exactly). (Hints: How are singular values of A related to eigenvalues of A^{t}A? What happens if you apply the 1-norm to an equation A^{t}Az = λz?)

Solution: The eigenvalues of A^{t}A are the squares of the singular values of A by Exercise 5.2.17. As in the Hint, let z be an eigenvector of A^{t}A with eigenvalue σ_{1}^{2}, and assume z has been normalized so ∥z∥_{1} = 1. Then ∥A^{t}Az∥_{1} = ∥σ_{1}^{2}z∥_{1} = σ_{1}^{2}. Hence

σ_{1}^{2} ≤ max_{∥x∥_{1}=1} ∥A^{t}Ax∥_{1} = ∥A^{t}A∥_{1} ≤ ∥A^{t}∥_{1}∥A∥_{1} = ∥A∥_{∞}∥A∥_{1},

(since the matrix 1-norm is the “maximum column sum” and the matrix ∞-norm is the “maximum row sum.”) Taking square roots gives

∥A∥_{2} = σ_{1} ≤ \sqrt{∥A∥_{1}∥A∥_{∞}}

which is what we wanted to show.

C) (10) Show using MATLAB that the inequalities in parts A and B are satisfied for the matrix B − I_{60}, where B is the 60 × 60 bucky matrix.

Solution: For part A, using MATLAB, we find the matrix B − I_{60} has rank 51, so

∥B − I_{60}∥_{F} = 15.4919 ≤ 25.8379 = \sqrt{51} · ∥B − I_{60}∥_{2}.

For part B,

∥B − I_{60}∥_{2} = 3.6180 ≤ 4.0 = \sqrt{∥B∥_{1}∥B∥_{∞}}.

II. More on condition numbers. Let n > m. Recall that if A ∈ M_{n×m}(R), the condition number κ_{2}(A) = σ_{1}/σ_{m} measures the susceptibility of the least-squares solution of Ax = b to round-off errors.

A) (10) Show that if an additional column y ∈ R^{n} is appended to A, to yield

\overline{A} = ( A y) ∈ M_{n×(m+1)}(R),

then

σ_{1} (\overline{A}) ≥ σ_{1} (A) and σ_{m+1} (\overline{A}) ≤ σ_{m} (A) .

What does this say about κ_{2}(\overline{A}) vs. κ_{2}(A)?

Solution: Let v_{1} be a singular vector of A corresponding to the largest singular value (that is, column 1 of the V matrix in an svd). We know Av_{1} = σ_{1}u_{1} and v_{1}, u_{1} are

unit vectors. Now consider the vector \overline{v} = ( \begin{matrix}v_{1} \\ 0 \end{matrix} ) ∈ R^{m+1}. Since v_{1} is a unit vector

in R^{m}, this is also a unit vector. We have

∥\overline{A}\overline{v}∥_{2} = ∥(A y) ( \begin{matrix}v_{1} \\ 0 \end{matrix} ) ∥_{2} = ∥Av_{1}∥_{2} = σ_{1}(A).

Hence

σ_{1}(\overline{A}) = ∥\overline{A}∥_{2} = max_{∥\overline{x}∥_{2}=1} ∥A\overline{x}∥_{2} ≥ σ_{1}(A),

where the maximum here is over all unit vectors in R^{m+1}. Similarly, if v_{m} is a singular vector of A corresponding to the smallest singular value (that is, column m of the V matrix in an svd). We know Av_{m} = σ_{m}u_{m} and v_{m}, u_{m}

are unit vectors. Now consider the vector \overline{v} = ( \begin{matrix}v_{m} \\ 0 \end{matrix} ) ∈ R^{m+1}. Since v_{m} is a unit

vector in R^{m}, this is also a unit vector. We have

∥\overline{A}\overline{v}∥_{2} = ∥(A y) ( \begin{matrix}v_{m} \\ 0 \end{matrix} ) ∥_{2} = ∥Av_{m}∥_{2} = σ_{m}(A).

Hence

σ_{m+1}(\overline{A}) = min_{∥\overline{x}∥_{2}=1} ∥A\overline{x}∥_{2} ≤ σ_{m}(A),

since the minimum here is over all unit vectors in R^{m+1}. (See Exercise 4.2.5 in the text – that is stated for square matrices, but the same is true in general.) Finally, combining the two inequalities shown here,

κ_{2}(\overline{A}) = \frac{σ_{1}(\overline{A})}{σ_{m+1}(\overline{A})} ≥ \frac{σ_{1}(A)}{σ_{m}(A)} = κ_{2}(A).

B) (10) Show that if an additional row, w^{t} for w ∈ R^{m}, is appended to A to yield

\overline{A} = ( \begin{matrix}A \\ w^{t} \end{matrix} ) ∈ M_{(n+1)×m}(R),

then

σ_{1} (\overline{A}) ≤ \sqrt{σ_{1}(A)^{2} + ∥w∥_{2}^{2}} and σ_{m} (\overline{A}) ≥ σ_{m} (A) .

What does this say about κ_{2}(\overline{A}) vs. κ_{2}(A)?

Solution: We have, taking maximum over vectors x ∈ R^{m},

σ_{1}(\overline{A}) = max_{∥x∥_{2}=1} ∥\overline{A}x∥_{2}

= max_{∥x∥_{2}=1} ∥ ( \begin{matrix}A \\ w^{t} \end{matrix} ) x∥_{2}

= max_{∥x∥_{2}=1} \sqrt{∥Ax∥_{2}^{2} + (w^{t}x)^{2}}

Now, by Cauchy-Schwarz, since x is a unit vector,

(w^{t}x)^{2} = ⟨w, x⟩^{2} ≤ ∥w∥_{2}^{2}∥x∥_{2}^{2} = ∥w∥_{2}^{2}.

Moreover, max ∥Ax∥_{2}^{2} = ∥A∥_{2}^{2} = σ_{1}(A)^{2}. Hence,

σ_{1}(\overline{A}) ≤ \sqrt{σ_{1}(A)^{2} + ∥w∥_{2}^{2}},

as claimed. Using the same formulas as above, for all unit vectors x ∈ R^{m},

∥\overline{A}x∥_{2} = \sqrt{∥Ax∥_{2}^{2} + ⟨w, x⟩^{2}} ≥ ∥Ax∥_{2}.

Hence the minimum over all x of ∥\overline{A}x∥_{2} must also be greater than or equal to the minimum over all x of ∥Ax∥_{2}. This shows σ_{m}(\overline{A}) ≥ σ_{m}(A).

Unlike the situation in part A, we cannot say anything about the relative sizes of κ_{2}(\overline{A}) and κ_{2}(A) here. Either one could be larger.

III. An image-processing application of the SVD. A large, detailed, image stored as a matrix of gray-scale pixel values can take a large amount of storage space. If the information contained in such an image can be compressed without losing too much image quality, that is a good thing for transmission and storage. One possible method for image compression is based on the the theoretical result on SVD’s in part A below. The later parts will show you how this works in practice.

A) (10) Recall from Theorem 4.1.12 (Exercise 4.1.13) in Watkins that if A = UΣV^{t} is an SVD of A, then if A has rank r,

(∗)  A = ∑_{j=1}^{r} σ_{j}u_{j}v_{j}^{t},

where u_{j} and v_{j} are the columns of the U and V matrices respectively, and σ_{1} ≥ σ_{2} ≥ ⋯ ≥ σ_{r} > 0 as usual. If we keep only the largest k singular values for some k < r, then the resulting matrix

(∗∗)  A_{k} = ∑_{j=1}^{k} σ_{j}u_{j}v_{j}^{t},

can be thought of as a rank k approximation to A (this makes especially good sense if the omitted singular values σ_{k+1}, …, σ_{r} are “small” compared to the others). Prove that:

(1) ∥A − A_{k}∥_{2} = σ_{k+1}, and

Solution: From (*) and (**) above,

A − A_{k} = ∑_{j=k+1}^{r} σ_{j}u_{j}v_{j}^{t}.

It follows that the singular values of A − A_{k} are σ_{k+1}, …, σ_{r} and the rest zero. Since σ_{k+1} ≥ σ_{j} for j > k + 1, this says ∥A − A_{k}∥_{2} = σ_{k+1}.

(2) if A′ is any other matrix of rank ≤ k, then

∥A − A′∥_{2} ≥ ∥A − A_{k}∥_{2}.

In words, A_{k} is the closest matrix to A among all matrices of rank k or less. Hint: For part (2), show that

N (A′) ∩ Span{v_{1}, …, v_{k+1}} ≠ {0}

and see what A − A′ does to a unit vector in the intersection.

Solution: If rank(A′) ≤ k, then by the fundamental equation

dim rank(A′) + dim N (A′) = m,

we see dim N (A′) ≥ m − k. Now V = N (A′) and W = Span{v_{1}, …, v_{k+1}} are two subspaces of R^{m} with dim(V) + dim(W) ≥ (m − k) + (k + 1) = m + 1. It follows by a general linear algebra fact in this situation that V ∩ W ≠ {0} (that is, there is some nonzero vector in the intersection of the two subspaces). Take any such z ∈ V ∩ W. We can normalize z to obtain a unit vector (divide z by ∥z∥_{2}). Then since z ∈ W, we have

z = c_{1}v_{1} + ⋯ + c_{k+1}v_{k+1}

and at least one c_{i} ≠ 0. But now, since z ∈ V = N (A′) as well,

(A − A′)z = Az − A′z = Az − 0 = Az.

Then, since v_{1}, …, v_{k+1} and u_{1}, …, u_{k+1} are orthogonal unit vectors, by the Pythagorean Theorem:

∥(A − A′)z∥_{2} = ∥Az∥_{2} = ∥A(c_{1}v_{1} + ⋯ + c_{k+1}v_{k+1})∥_{2}
= ∥c_{1}σ_{1}u_{1} + ⋯ + c_{k+1}σ_{k+1}u_{k+1})∥_{2}
= \sqrt{(c_{1}σ_{1})^{2} + ⋯ + (c_{k+1}σ_{k+1})^{2}}
≥ σ_{k+1}\sqrt{c_{1}^{2} + ⋯ + c_{k+1}^{2}}
= σ_{k+1}∥z∥_{2}
= σ_{k+1}.

This implies ∥A − A′∥_{2} ≥ σ_{k+1}, which, when combined with part (1), is what we wanted to show.

B) (5) One measure of the size of an image A is the total number of real numbers needed to write the vectors u_{j}, v_{j} and the σ_{j} in (*) or (**). If A is 200 × 320 and has rank 200, what is the

compression ratio = \frac{size using (∗∗)}{size using (∗)}

achieved if we replace the original expression (*) for A with (**), using k = 5, 10, 20, 25?

Solution: With k = 5:

\frac{5 + 5 × 200 + 5 × 320}{200 + 200 × 200 + 200 × 320} = .025

With k = 10:

\frac{10 + 10 × 200 + 10 × 320}{200 + 200 × 200 + 200 × 320} = .05

With k = 20:

\frac{20 + 20 × 200 + 20 × 320}{200 + 200 × 200 + 200 × 320} = .1

With k = 25:

\frac{25 + 25 × 200 + 25 × 320}{200 + 200 × 200 + 200 × 320} = .125

Comment: A fairer comparison might actually be to take the numerator here over the size of the matrix form of A. For instance with k = 25:

\frac{25 + 25 × 200 + 25 × 320}{200 × 320} = .2035

By that measure, the “compressed” form with k = 25 is about 1/5 the size of the original image in the matrix format.

C) (10) Using MATLAB, test out this compression scheme using the image file clown.mat:

load clown.mat;
colormap('gray');

This will store the image file as a 200 × 320 full matrix called X. You can display the full image using

image(X)

Now, compute the SVD of X calling the factors U, S, V. To compute the “compressed” matrices A_{k} for various k, you can use commands like this:

U(:,1:k)*S(1:k,1:k)*V(:,1:k)'

(you supply the values of k). Note: the MATLAB syntax A(:, a : b) means: form the submatrix of A taking all rows and columns a through b.) For k = 5, 10, 20, 25, compute A_{k}, display the resulting images, and comment on how well they represent the full image (k = 200). As a more precise measure of image quality, also compute ∥A − A_{k}∥_{2} for each of these k values.

Solution: You should have seen that the image quality (judging by eye) steadily improved with k = 5, 10, 20, 25, until the image with k = 25 was hardly distinguishable from the original. The precise measures of image quality are (all ×10^{3}):

\begin{aligned}∥A − A_{5}∥_{2} & = σ_{6} = 1.0699 \\ ∥A − A_{10}∥_{2} & = σ_{11} = .6250 \\ ∥A − A_{20}∥_{2} & = σ_{21} = .3288 \\ ∥A − A_{25}∥_{2} & = σ_{26} = .2610 \end{aligned}

IV. More on iterative methods. Recall that the Jacobi and Gauss-Seidel iterative methods for square systems Ax = b can be derived by splitting the coefficient matrix A as a sum. The general idea would be to write A = M + N for some square matrices M, N with M invertible.

A) (5) Show that however this is done, the resulting iteration can be written as a correction based on the residual r^{(k)} = b − Ax^{(k)}:

x^{(k+1)} = x^{(k)} + M^{−1}r^{(k)}.

Solution: The rearrangement to fixed-point form based on the splitting A = M + N is found by starting from Ax = (M + N)x = b:

x = −M^{−1}Nx + M^{−1}b

This leads to the iteration formula

x^{(k+1)} = −M^{−1}Nx^{(k)} + M^{−1}b.

Now, if we rewrite N as A − M, this becomes

x^{(k+1)} = −M^{−1}(A − M)x^{(k)} + M^{−1}b = x^{(k)} + M^{−1}(b − Ax^{(k)}) = x^{(k)} + M^{−1}r^{(k)},

which is what we wanted to show.

B) (5) For the remainder of this problem, assume that all eigenvalues of A are real and non-negative. The method obtained with M = \frac{1}{ω} I for some ω > 0 and N = A − M is called Richardson’s method. Richardson’s method iteration in the fixed point form is

x^{(k+1)} = (I − ωA)x^{(k)} + ωb.

Show that Richardson iteration converges only for ω < \frac{2}{λ_{max}}, where λ_{max} is the largest eigenvalue of A.

Solution: This is similar to one problem from Problem Set/Lab 11. The Richardson iteration converges if and only if the spectral radius of I − ωA is < 1. The spectrum of this matrix is the set

σ(I − ωA) = {1 − ωλ : λ ∈ σ(A)}.

Since all λ ≥ 0 and ω > 0,

σ(I − ωA) ⊂ [1 − ωλ_{max}, 1 − ωλ_{min}] ⊂ (−∞, 1) ⊂ R

for all ω. The spectral radius is < 1 only if

1 − ωλ_{max} > −1 ⟺ ω < \frac{2}{λ_{max}} .

C) (10) Show that the omega that minimizes the spectral radius of the “Richardson

G-matrix” I − ωA is

ω_{opt} = \frac{2}{λ_{max} + λ_{min}} ,

where “opt” stands for “optimal – explain why this would be the best value of ω to use.

Solution: As a function of ω, the spectral radius of I − ωA is given by the maximum of |1 − ωλ_{max}| and |1 − ωλ_{min}|. Which one is larger depends on whether 1 − ωλ_{max} is closer to -1 or 1 − ωλ_{min} is closer to 1.

\begin{aligned}ρ(I − ωA) & = max\{|1 − ωλ_{max}|, |1 − ωλ_{min}|\} \\ & = \{ \begin{matrix}|1 − ωλ_{max}| & if ωλ_{min} < 2 − ωλ_{max} \\ |1 − ωλ_{min}| & if ωλ_{min} > 2 − ωλ_{max} \end{matrix} \end{aligned}

Plotting the two parts of this function separately on the interval (0, 2/λ_{max}), we see that both start at 1 with ω = 0. The first is a vee-shaped graph that goes down to zero at ω = 1/λ_{max} and comes back up to 1 at ω = 2/λ_{max}. The second is (part of) a second vee-shaped curve sloping down more gradually from 1 at ω = 0, possibly hitting the ω axis, and coming back up to |1 − 2λ_{min}/λ_{max}| at ω = 2/λ_{max}. Note: this could happen before you hit the “vee,” depending on the exact values of λ_{max} and λ_{min}. The smallest max occurs when the two lines intersect, which is when

ωλ_{min} = 2 − ωλ_{max} ⟺ ω = \frac{2}{λ_{max} + λ_{min}} .

This value is optimal in the sense that the spectral radius is minimized. The smaller the largest eigenvalue is (in absolute value), the more rapid the convergence will be.

D) (5) Refer to the system from Example 7.2.3 in the text. Using MATLAB, determine ω_{opt} for Richardson on this system, and determine the number of Richardson iterations needed to yield a solution that is accurate to 8 decimal places using ω = .17, ω_{opt}, and ω = .1.

Solution: We have λ_{max} = 11.2561, and λ_{min} = 3.3182, so

ω_{opt} = \frac{2}{11.2561 + 3.3182} = .1372.

The exact solution is (4, 3, 2, 1)^{t}. Iterations needed for 8 decimal place accuracy:

ω = .17 about 224
ω = ω_{opt} = .1372 about 34
ω = .1 about 48.