Math 550A MATLAB Assignment #4 5
P4 = downsamp(2); D2 = fftdiag(2); F2 = fft(eye(2));
F4H = [F2, D2*F2; F2, -D2*F2]*P4
Compare F4H with the matrix fft(eye(4))’. Are they the same?
Now repeat the process to generate the 8 × 8 matrix F
H
8
. First generate P8 = downsamp(4) and D4 =
fftdiag(4). Use these, the matrix F4H you just generated, and equation (29) in the supplementary notes
to obtain a matrix F8H. Compare your F8H with the Matlab-generated matrix F
H
8
by calculating norm(F8H
- fft(eye(8))’). This should be of size 10
−15
, thus zero to the limits of machine computation.
(b) Speed Comparison: Generate the 2
12
× 2
12
Fourier matrix and a random vector x ∈ R
4096
(be sure to include the ; at the end of the line so that the matrix will not be displayed or
written in your diary file):
F = fft(eye(4096))’; x = rand(4096,1);
Now calculate the computation time to obtain Fourier transform of x in two ways: first by direct matrix
multiplication and then by the fast Fourier transform (be sure to type the semicolons as indicated):
tic; F’*x; matrixtime = toc
tic, fft(x); ffttime = toc
speedup = floor(matrixtime/ffttime)
Now compare your observed speedup ratio using the FFT with the theoretical prediction. If N = 2
k
, then
computing the matrix × vector F’*x requires N
2
= 2
2k
scalar multiplications. Computing fft(x) needs at
most k2
k−1
scalar multiplications by Equation (30) in the supplementary class notes (also see Strang pp.
193-4). The theoretical speedup ratio for the FFT over the plain (matrix × vector) multiplication, in terms
of the scalar multiplications required, is thus
2
2k
/[k2
k−1
] = 2
k+1
/k . (?)
A similar speed advantage holds for the scalar additions required. Compute the theoretical speed advantage
(?) for the value k = 12 that you have used, and compare it with the observed value speedup. Scalar
multiplications require much more computing time than scalar additions, so it is the reduction in the number
of multiplications that is crucial for the speedup given by the FFT.
Question 5. The QR Eigenvalue Algorithm
“The QR algorithm is one of the most important, widely used, and successful tools we have in technical
computation. Several variants of it are in the mathematical core of Matlab. They compute the eigenvalues
of real symmetric matrices, real nonsymmetric matrices, and pairs of complex matrices, and the singular
values of general matrices.” (from Numerical Computing with Matlab, by Cleve Moler, the original author
of Matlab).
(a) Single-shift QR algorithm : This algorithm is described on page 364 of Strang (equation (6)).
To illustrate this algorithm, use the integer matrix A = gallery(3) from the Matlab files (this matrix has
eigenvalues 1, 2, and 3). Set n = size(A,1), I = eye(n,n) and then iterate the following Matlab line
using the up-arrow key:
s = A(n,n); [Q, R] = qr(A - s*I); A = R*Q + s*I
(i) Show by a hand calculation that the new matrix A produced by this line is orthogonally
similar to the old matrix A.
Hence the new A has the same eigenvalues as the old A (although the eigenvectors will be different). Since
the similarity is by an orthogonal transformation, the procedure is numerically stable.
Continue forming a new A from the old A until the last row of A becomes [0 0 2.000] (this should
require about ten iterations).
(b) Deflation: Now replace the matrix A by the 2 × 2 matrix A = A(1:2,1:2). Update the values of
n and I and repeat the QR iteration from part (a) until this 2 × 2 matrix is upper triangular (this should
require only two or three iterations).