English homework help. MATH466/462 Project 4. Due in class on Wed, Mar 18, 2020

Instruction: your project report should include necessary mathematical

justification, description, and details of your algorithms/conclusions.

Your MATLAB codes and generated outputs may be attached in the end of the

report. Make sure you addressed all the questions in each problem. Both

the report and codes will be graded. Please submit a printed hard-copy.

Problem A (20 pts): Circulant Preconditioners for Toeplitz Systems

Toeplitz matrix arises in many different applications. Its matrix-vector product can be computed efficiently via fast

Fourier transform (FFT). Many circulant preconditioners have been proposed for solving Toeplitz systems. In this

project, we will compare several circulant preconditioners for solving various Toeplitz systems.

Task 1: Study Toeplitz matrix, Circulant matrix, and the use of FFT (read our Textbook Chap 10.1).

Read Chapters 1 and 3: https://ee.stanford.edu/~gray/toeplitz.pdf

Understand more on FFT: https://arxiv.org/pdf/1805.05533v2.pdf

Task 2: play and understand the following codes:

For any Circulant matrix ?, the matrix-vector products ? = ?

−1

? can also be computed via FFT:

1 n=10; C=gallery gallery gallery(‘circul’,(1:n));%Construct a Circulant matrix using (1:n) as 1st row

2 b=ones(n,1); x2=C\b; %Direct solve: O(n^3) operations

3 ev=fft(C(:,1));%the eigenvalues of C by FFT of its first column

4 x1= ifft(fft(b)./ev); %Solve C\b using FFT: O(n log n) operations

5 norm(x1-x2,inf) %should be zero

By embedding a Toeplitz matrix ? into a Circulant matrix, the product ? ? can also be computed via FFT:

1 n=10;t=(n:-1:1); T=toeplitz toeplitz toeplitz(t,t’); %construct a full symmetric Toeplitz matrix: T’=T

2 v=rand(n,1);y2=T*v;%compute y2=T*v using direct multiplication: O(n^2) operations

3 gev = fft([t 0 t(n:-1:2)].’);%the eigenvalues of the embeding larger Circulant matrix

4 y = ifft(fft([v;zeros(n,1)]).*gev);%compute y1=T*v using FFT: only O(n log n) operations

5 y1 = y(1:n); %take the first half of the long vector

6 norm(y1-y2,inf) %should be close to zero

Task 3: understand the construction of 3 circulant preconditioners.

Let ?? be an ?-by-? Toeplitz matrix with ?? (?, ?) = ??−?

, where {?? }

?−1

?=1−?

are given diagonals. We can define at least

3 different circulant preconditioners as follows:

1. Strang’s Preconditioner: Strang’s preconditioner ?? with ?? (?, ?) = ??−?

is defined to be the circulant matrix obtained by copying the central diagonals of ?? and bringing them around to complete the circulant

requirement. Assume ? = 2? is even, the diagonals ?? of ?? are given by

?? =

?? if 0 ≤ ? ≤ ? − 1

0 if ? = ?

??−? if ? < ? ≤ ? − 1

?−? if (1 − ?) ≤ ? < 0

.

2. T. Chan’s Preconditioner: T. Chan’s preconditioner ?? with ?? (?, ?) = ??−?

is defined through minimizing

the difference between?? and?? over all circulat matrices. The diagonals?? of?? are given by (taking ?−? = 0)

?? =

{

(?−?)?? +???−?

?

if 0 ≤ ? ≤ ? − 1

??+? if (1 − ?) ≤ ? < 0

.

3. R. Chan’s Preconditioner: R. Chan’s preconditioner ?? with ?? (?, ?) = ??−?

is defined to make uses of all

the entries of ??. The diagonals ?? of ?? are given by (taking ?−? = 0)

?? =

{

?? + ??−? if 0 ≤ ? ≤ ? − 1

?−? if (1 − ?) ≤ ? < 0

.

1

(1) Use our PCG (mypcgfun.m) to solve the SPD Toeplitz systems ??? = ? with ? = ????(?, 1) and

?? (?, ?) =

{

?

2

/3 if ? == ?

2(−1)

?−?

(?−?)

2 if ? ≠ ?

.

Test with no preconditioner and the above circulant preconditioners, and compare their iteration numbers

and CPU times for different dimensions n=1e3*(1:5) (set max 10000 iterations and tolerance 10−7

).

You can start with construction of all related full matrices for smaller ?, but eventually all matrix-vector

products should be replaced by only FFT based codes for better efficiency and lower memory costs.

(2) Circulant preconditioners are attrative since they can be solved efficiently (?(? ln?) operations) via FFT. As we

aready know, a tridiagonal matrix can also be solved very efficiently (?(?) operations) via Thomas algorithm.

Run your codes again with the following tridiagonal preconditioner, which works well for the given ??:

1 Pn=gallery gallery gallery(‘tridiag’,n,-1,2,-1);

For benchmarking your own codes, below are the iteration numbers based on my implementation:

1 n None Strang T. Chan R. Chan Tridiag

2 1000 746 8 28 7 14

3 2000 1522 8 35 7 14

4 3000 2301 8 41 7 14

5 4000 3081 8 47 7 14

6 5000 3862 8 50 7 14

(3) Solve the same system using the Gaussian elimination method (GEsolver.m), what you observe in terms

of CPU times growth, in comparison with the above PCG solvers?