MATH3204 Numerical Linear Algebra & Optimisation (S2-2022): Assignment 02
MATH3204/7234 (S2-2022): Assignment 02
1. [5 marks each]
(a) Show that if A∈C and λ is an eigen value of A, then λ is an eigen value of A*.
(b) Let A ∈ C . For an arbitrary matrix induced norm show that
ρ(A) ≤ ||Ak ||1/k, k = 1,2,...,
where ρ(A) is the spectral radius of A.
(c) Consider A, B ∈ C . Prove that if one of A or B is non-singular, then AB and BA are similar.
2. [10 marks each] SupposeA∈Cm×n, B∈Cn×m,andm≤n.
1. (a) We know AB and BA are not in general similar. However, we can still say some- thing about their spectrum. Show that the n eigenvalues of BA are the m eige- navlues of AB together with n − m zeros.
2. (b) Suppose A ∈ C is invertible and x, y ∈ C . Use the first part of this question and show that
det(A + xyT) = det(A)(1 + yTA−1x).
3. [5 marks]
Consider the least squares problem
min 1 ∥Ax − b∥2 ,
x∈Rn 2
where A ∈ Rm×n is rank deficient. Show that among all possible solutions to this problem, A†b has the smallest Euclidean norm.
4. [5 marks each]
In this question, you will study one of the reasons why forming the matrix-inverse explicitly is so much frowned upon.
(a) For n = 100,200,500,1000,2000,5000,10000, consider a series of linear systems Anxn = bn where An = randn(n, n), bn = ones(n, 1) if you are coding in Matlab or An = np.random.randn(n, n), bn = np.ones((n, 1)) in Python.
(b) Does your observation in comparing the time for these two implementations roughly match that predicted by the theory in terms of “flops”? Provide some explanation for your answer.
On paper, these functions are implementing the same thing, namely x = A−1b , but on computer, they perform differently! This is one of the places where numerical linear algebra differs from linear algebra.
(c) Now generate a similar plot and make similar comparisons as above. However, this time construct An using sprandn(n,n,0.1) or sp.sparse.random(n,n,0.1, format=’csc’) if you are using Matlab or Python, respectively. These commands construct random but sparse matrices. Feel free to check the relevant documenta- tion and see what each argument of these functions is used for. The third argument controls the density of the matrix, i.e., what percentage of the entries of this matrix consists of non-zeros values.
(d) Does your observation in comparing the time for these two implementations roughly match that predicted by the theory in terms of “flops”? Provide some explanation for your answer.
5. [5 marks each] For solving Ax = b, consider the stationary Richardson iteration
xk+1 =xk +α(b−Axk).
For this question, we assume that all eigenvalues of A are real.
1. (a) Find the iteration matrix.
2. (b) Show that if λmin < 0 < λmax, the method will always be divergent for some initial iterate.
3. (c) So assuming that all eigenvalues are positive, find the spectral radius of the iteration matrix in terms of α, λmin and λmax ?
4. (d) Find αmin < αmax such that the method converges for any αmin < α < αmax. In other words, find the largest range for α that ensures the method always stays convergent.
5. (e) What is the optimal value for α? In other words find αopt that minimizes the spectral radius of the iteration matrix.
6. (f) Suppose A ≻ 0 . Using this optimal step-size, obtain the optimal convergence rate in terms of the matrix condition number.
7. (g) Implement the stationary Richardson iteration and apply it to the matrix gener- ated using the following code. For this assignment, we use beta = gamma = 0, N = 50, x0 = 0, and a right hand-side vector b containing all ones. Compare the performance using three step-sizes: αopt,0.5αopt, and 1.1αopt. What do you observe? Terminate the algorithms if the Euclidean norm of the residual vector ∥rk∥2 = ∥b − Axk∥2 ≤ 10−12 or if a maximum number of iteration of 20, 000 is reached.
function A=lcd(beta ,gamma ,N)
%
% Simple code to generate test matrices for this question and
maybe some later assignments % Usage:
% to generate symmetric positive definite, call with beta= gamma =0
% to generate nonsymmetric call, e.g., with beta=gamma=1/2 %
% Note: output matrix is of size N^2-by-N^2
ee=ones(N,1);
a=4; b=-1-gamma; c=-1-beta; d=-1+beta; e=-1+gamma;
t1=spdiags([c*ee,a*ee,d*ee],-1:1,N,N);
t2=spdiags([b*ee,zeros(N,1),e*ee],-1:1,N,N);
A=kron(speye(N),t1)+kron(t2,speye(N));
end
For Matlab use the above code snippet.
import numpy as np
import scipy as sp
def lcd(beta ,gamma ,N):
# Simple code to generate test matrices for this question and maybe some later assignments
# Usage:
# to generate symmetric positive definite, call with beta =gamma=0
# to generate nonsymmetric call, e.g., with beta=gamma =1/2
#
# Note: output matrix is of size N^2-by-N^2
ee = np.ones((1, N))
a, b, c, d, e = 4, -1 - gamma, -1 - beta, -1 + beta, -1 + gamma
t1 = sp.sparse.spdiags(np.vstack([c * ee, a * ee, d * ee]), np.arange(-1,2), N, N)
t2 = sp.sparse.spdiags(np.vstack([b * ee, np.zeros((1, N)), e * ee]), np.arange(-1,2), N, N)
A = sp.sparse.kron(sp.sparse.eye(N,N), t1) + sp.sparse.kron( t2, sp.sparse.eye(N,N))
For Python use the above code snippet.
(h) Implement the accelerated variant of the Richardson iteration and, on the same test problem and starting from the same x0 as above, compare its performance with the stationary version. For both methods, use the step-size α = 2/(λmin + λmax).
100 marks in total
Note:
- This assignment counts for 15% of the total mark for the course.
- Although not mandatory, if you could type up your work, e.g., LaTex, it would be
greatly appreciated.
- Show all your work and attach your code and all the plots (if there is a programming question).
- Combine your solutions, all the additional files such as your code and numerical results, all in one single PDF file.
- Please submit your single PDF file on Blackboard.
