matlab code

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