Numerical Methods 2021 — Midterm 2
Solutions
Problem 1 (20 pts). Consider an upper triangular matrix A stored in a compressed row format (CSR)
(a) (5 pts) Let
A =
5 2 8 0 4 1 0 0 3
Show how A is stored in compressed row format in the following table.
Here we consider array index beginning with 1.
For column indexes with the same row index, you should store them in an ascending order in the table.
array index 1 2 3 4 5 6 a
acol ind arow ptr
(b) (10 pts) Assume A is stored in the same way as in (a). Give a code to solve any upper triangular linear equation
Ax = b, where we assume Aii 6= 0.
(c) (5 pts) Let
A =
5 2 8 0 4 1 0 0 3
and b =
68 11 21
Run your code to solve the system and output your answer.
Solution.
(a) The CSR matrix format of A is
array index 1 2 3 4 5 6
a 5 2 8 4 1 3
acol ind 1 2 3 2 3 3 arow ptr 1 4 6 7
(b) function x = upper_linear_solver(A, b)
n = size(A, 2); % A is m, n dim array [a, acol_ind, arow_ptr] = getCSR(A);
for i = n:-1:1
Aijxj = 0; % sum of Aij * xj, where j = i+1, ..., n Aii = a(arow_ptr(i));
for k = (arow_ptr(i)+1):(arow_ptr(i+1)-1) j = acol_ind(k); % column c, value v Aijxj = Aijxj + a(k) * x(j);
end
x(i) = (b(i) - Aijxj) / Aii;
end end (c) iter 1
Aijxj = 0
x[3] = b[3] - Aijxj / 3 = 7 iter 2
Aijxj = 0
Aijxj = Aijxj + a[5] * x[acol_ind[5]]
= Aijxj + (a[5] * x[3])
= Aijxj + (1 * 7)
= 7
x[2] = b[2] - Aijxj / 4 = 1 iter 3
Aijxj = 0
Aijxj = Aijxj + a[2] * x[acol_ind[2]]
= Aijxj + (a[2] * x[2])
= Aijxj + (2 * 1)
= 2
Aijxj = Aijxj + a[3] * x[acol_ind[3]]
= Aijxj + (a[3] * x[3])
= Aijxj + (8 * 7)
= 58
x[1] = b[1] - Aijxj / 5 = 2 ans =
2 1 7
Problem 2 (30 pts). Consider the following linear equation, Ax = b,
where A is an n-by-n symmetric positive-definite matrix. Assume A is a sparse matrix stored in com- pressed column (CSC) format. For row indexes with the same column index, we require that the CSC format stores them in an ascending order.
(a) (10 pts) In the Gauss-Seidel method, in each iteration we update the i-th element of x, xi, by a value δ
xi ← xi+ δ.
Suppose we already cached a vector
c = b − Ax,
How to represent the value δ by using values in c and diagonal values of A?
Suppose we would like to maintain the vector c by using the new x. How to update c by using δ and components in A?
(b) (15 pts) Give the code to do the Gauss-Seidel method by using the technique in the prob. (a).
We require that it can take any vector as the first iterate.
Note that you may need to extract diagonal elements to a dense vector first.
Also, you may need to initialize an dense vector c to adapt the technique derived in prob. (a).
(c) (5 pts) Let
A =
4 1 1 1 4 0 1 0 3
and b =
18 25 25
Show how A is stored in compressed column format (CSC) in the following table.
Here we consider array index beginning with 1.
1 2 3 4 5 6 7 a
arow ind acol ptr Given initial solution
x0 =
0 14
0
Run your code to solve the system. You need to details such as δ, x and the vector c at each step.
Solution.
(a) The new value ¯xi by using the Gauss-Seidel method is
¯
xi = bi−P
j6=iAijxj
Aii . If we minus both sides with the original value xi, we then have
δ = ¯xi− xi = bi−P
j6=iAijxj
Aii − xi = bi−P
jAijxj Aii = ci
Aii. For c, we update it by
c = b − A(x + δei) = c − δ
A1i
... Ani
, where ei = [0, · · · , 1, · · · , 0]T is a standard unit vector.
(b) function y = Gauss_Seidel(A, x, b) n = size(A, 2);
[a, arow_ind, acol_ptr] = getCSC(A);
% Initialize b - Ax and Adiag for i = 1:n
c(i) = b(i);
for j = acol_ptr(i):(acol_ptr(i+1)-1) r = arow_ind(j); v = a(j);
c(r) = c(r) - v * x(i);
if i == r
Adiag(i) = v;
end end end
% Run Gauss_Seidel until b - Ax is zero vector while sum(c) ˜= 0
for i = 1:n
delta = c(i) / Adiag(i);
if delta ˜= 0
x(i) = x(i) + delta;
for j = acol_ptr(i):(acol_ptr(i+1)-1) r = arow_ind(j); v = a(j);
c(r) = c(r) - v * delta;
end end end end y = x;
end
(c) The CSC format of A is
1 2 3 4 5 6 7
a 4 1 1 1 4 1 3
arow ind 1 2 3 1 2 1 3 acol ptr 1 4 6 8
iter 1
update x[1]
delta = 1
x = [1 14 0]
c = [0 -32 24]
update x[2]
delta = -8
x = [1 6 0]
c = [8 0 24]
update x[3]
delta = 8
x = [1 6 8]
c = [0 0 0]
y = 1 6 8
Problem 3 (35 pts). Consider a linear system Ax = b:
2 0 2 0 1 0 2 0 3
x =
−1
−1 1
(a) (5 pts) Is the matrix A symmetric positive definite? If True, please prove it. Otherwise, give a counterexample to show that is False.
(b) (5 pts) Take x0 =0 0 0T. Do two CG iterations and show what x1 and x2 are.
(c) (5 pts) Is x2 a solution or not?
(d) (10 pts) Let r1 and p1 be vectors from the aforementioned procedures. Solve the following opti- mization problem with the variable p:
minp kp − r1k22
s.t. p ∈ span{Ap1}⊥
(1)
How is the solution of this problem connected to p2 obtained in the CG procedure?
(e) (5 pts) If we denote the solution of (1) as ¯p2 and calculate ¯α2, x2 by
¯
α2 = p¯T2r1
¯
pT2A¯p2, x2 = x1+ ¯α2p¯2. What are ¯α2, x2 ?
(f) (5 pts) In slides, we have a theorem by writing
A = I + B
and checking the relationship between rank(B) and the number of CG steps. From what you have observed in (a)-(d), can you say that rank(B) must be greater or equal to some value? Then check B to confirm the result.
One may indeed make mistakes in doing the calculation. However, if you understand the concept of CG, you should be able to easily validate your results.
Solution.
(a) Since
A = AT, A is symmetric. Given a vector v = [v1, v2, v3]T,
vTAv = 2v21+ 4v1v3+ 3v32+ v22 = 2(v1+ v3)2+ v23+ v22 That implies
vTAv = 0 iff v = 0.
Therefore, A is a symmetric positive definite matrix.
(b)
x0 =
0 0 0
, r0 = b =
−1
−1 1
,
⇒p1 =
−1
−1 1
, α1 = rT0r0 pT1Ap1 = 3
2, x1 = x0 + α1p1 =
−3/2
−3/2 3/2
, r1 = r0− α1Ap1 =
−1 1/2
−1/2
⇒β2 = rT1r1 rT0r0 = 1
2, p2 = r1+ β2p1 =
−3/2 0 0
, α2 = rT1r1 pT2Ap2 = 1
3, x2 = x1+ α2p2 =
−2
−3/2 3/2
,
(c) Since
r2 = r1 − α2Ap2 =
0 1/2 1/2
is not zero, x2 is not the solution.
(d) We can derive that
span{Ap1} =
0
−t t
, t ∈ R, so span{Ap1}⊥ is equivalent to
−p2 + p3 = 0 Therefore, (1) is equivalent to
minp
p1+ 1 p2− 1/2 p3+ 1/2
2
2
≡ min
p
p1+ 1 p2− 1/2 p2+ 1/2
2
2
≡ min
p (p1+ 1)2 + (p2− 1/2)2+ (p2+ 1/2)2 s.t. − p2+ p3 = 0
and
(p1+ 1)2+ (p2− 1/2)2+ (p2+ 1/2)2 = (p1+ 1)2+ 2p22+ 1/2.
Thereby, the solution is −1 0 0T, and we can see that p is parallel to p2.
(e)
¯ p2 =
−1 0 0
, ¯α2 = p¯T2r1 p¯T2A¯p2 = 1
2, x2 = x1+ ¯α2p¯2 =
−3/2
−3/2 3/2
+ 1 2
−1 0 0
=
−2
−3/2 3/2
. (f) From the theorem on page 6 of Conjugate gradient methods part4,
the rank of B is at least two (2)
because x2 is not a solution. Let us check on the rank of B. The matrix B can be found by
A =
2 0 2 0 1 0 2 0 3
= I +
1 0 2 0 0 0 2 0 2
= I + B
The rank of B is two since there are only two linear independent rows in B. Thus, the result (2) holds.
Problem 4 (15 pts). In ours slides, we use
f (x) = cos(x)
as an example to run the Newton method. Now we would like to check if it satisfies the theorem proved in slides and have the quadratic convergence.
(a) (5 pts) Is f (x) Lipschitz continuous? If so, prove it and give an α.
(b) (5 pts) For this particular f (x), can you directly apply Taylor expansion to prove Lemma 3? (Hint:
Check the Taylor theorem.) (c) (5 pts) Now consider
x∗ = π 2 and
{xk} → x∗.
Find ¯β and δ such that the second result of the theorem holds. That is,
|xk+1− x∗| ≤ δ|xk− x∗|2, ∀k ≥ L.
Solution.
(a) We know that
f0(x) = − sin(x).
By Mean Value Theorem, for any x and y, there exists t between x and y such that sin(x) − sin(y) = (x − y) sin0(t).
We have
| sin0(t)| = | cos(t)| ≤ 1, so
| sin(x) − sin(y)| = |(x − y) sin0(t)| = |x − y|| cos(t)| ≤ |x − y|, ∀x, y.
Thus we can choose
α = 1.
(b) Taylor Theorem tells us that the Taylor expansion of f (y) at x:
cos(y) = cos(x) − sin(x)(y − x) − 1
2cos(t)(y − x)2, where t is between x and y. Therefore, the error term is
e(x, y) = −1
2cos(t)(y − x)2, and
|e(x, y)| = |1
2cos(t)(y − x)2|
= 1
2| cos(t)||(y − x)2|
≤ 1
2|y − x|2 (c) Follow the slides, we have
β = |f¯ 0(x∗)−1| =
1 sin(π/2)
= 1 and
δ = α ¯β = 1.
