MATH466/46

2

Project 4. Due in class on Wed, Mar

1

8, 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
Understand more on FFT: https://arxiv.org/pdf/1805.05533v2
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=gallerygallerygallery(‘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=fftfftfft(C(:,1));%the eigenvalues of C by FFT of its first column
4 x1= ifftifftifft(fftfftfft(b)./ev); %Solve C\b using FFT: O(n log n) operations
5 normnormnorm(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=toeplitztoeplitztoeplitz(t,t’); %construct a full symmetric Toeplitz matrix: T’=T
2 v=randrandrand(n,1);y2=T*v;%compute y2=T*v using direct multiplication: O(n^2) operations
3 gev = fftfftfft([t 0 t(n:-1:2)].’);%the eigenvalues of the embeding larger Circulant matrix
4 y = ifftifftifft(fftfftfft([v;zeroszeroszeros(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 normnormnorm(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 ma-
trix 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

https://ee.stanford.edu/~gray/toeplitz

https://arxiv.org/pdf/1805.05533v2

(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=gallerygallerygallery(‘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?

2