18.330 Problem set 8

Exercise 1 : Eigenvalue solvers for a special matrix

In this problem we will consider a very special symmetric matrix. Recall that the second-order finite difference scheme for a function ( ) is given by

″

=

+1 − 2 + −1

ℎ0 .

= ( 0

+ )

ℎ =

2

−

where

and

Define the vector f such that f = . Then the corresponding derivative vector using finite differences can be written as f″ = f . Some care must be taken about the boundary conditions; for simplicity we will assume that 0 = = 0. Under these assumptions the matrix is a tridiagonal matrix:

′

−2

1

1

1′

2′

1−21

2

3

1

−2 1

3

′ ⋮

=

⋱ ⋱ ⋱

1

⋮

′

1

−2

−2

−2

1

−1

−2 −1

You can construct the matrix in Julia using the diagm function. For simplicity keep everything dense in this problem, although there clearly a lot of structure that could be exploited.

1. Consider the matrix . Using the ansatz (i.e. hypothesis) v = sin( / ), where 1 ≤ ≤ − 1, show that v is an eigenvector of . What is the corresponding eigenvalue? Remember that the eigenvalues should be real!

2. In class we discussed several ways to find eigenvalues of matrices. The simplest algorithm for finding eigenvalues is the power method. Supppse that a symmetric matrix has the eigendecomposition = Λ −1. Recall that for a symmetric matrix the eigenvectors associated with distinct eigenvalues are orthogonal. Starting with a random intitial vector x which can be written as a sum over the eigenvectors,

x = ∑ v

1

show that

= ∑

x v

where we order the eigenvalues by their magnitude, | 1| > | 2| > ⋯ > | | and hence

x = 1 1 (v1 + ( 2 v2)) .

1

This shows why the power method converges to the leading eigenvector.

̃

3. The power method gives an approximation for the leading eigenvector v1.

Since it is only an approximation it is not necessarily true that v ̃ =

1

̃

1v1 exactly.

Instead we will say that the best approximation to 1 is given by the that satifies

1 =

min

v

v

v

v

2

||

1̃−

1̃||2 = ∑ (∑ (

1̃)−(

1̃)) .

v

By differentiating this

expression with respect to

, show that

v

v

v

v

v

v⊤ v

1̃⋅ ( 1̃)

1̃ 1̃

1 ≈

=

⊤

.

1̃⋅

1̃

1̃

1̃

This is called the Rayleigh quotient.

4. Implement the power method power_method(A, N, x) for a symmetric matrix which iterates for a fixed number of iterations on the intial

vector x and returns ṽand approximated by the Rayleigh quotient

1 1

above.

5. Run this on the matrix 10 which is of size (9 × 9). Use the true largest magnitude vector from part 1. Plot the relative error between your result and the true value as a function of to get a smooth plot; use the same initial vector each time.

Remember to normalize the vector at each iteration to avoid overflow! This initial random vector should be complex. (Extra credit: show that the relationship is what you would expect analytically!)

2

6. A more advanced method that we discussed to find all the eigenvalues was the QR algorithm. We define (0) = and then

() ()= ()

(+1)= () ()

This will normally converge to an upper-triangular matrix. How does this work? We call two matrices and similar if = ⊤ where is orthogonal. Show that if has eigenvalue pairs ( , v ) then has eigenpairs ( , v ).

7. Show that in the QR algorithm ( +1) is similar to ( ) and hence . Therefore if the algorithm converges to an upper-triangular matrix we can read off the eigenvalues from the diagonal and the eigenvectors will be given by the columns of (1) (2) (3) ⋯ ( ).

8. Implement the QR algorithm in a function qr_algorithm(A, N) for a ma-trix with iterations. It should return the resulting matrix ( ) and the eigenvector matrix. Run the algorithm on a random symmetric matrix of size 10 × 10 for 100 iterations. How close to upper-triangular is the resulting matrix?

9. Run the QR algorithm on the matrix 10 of size (9 × 9) for a suitable number of iterations to get convergence. Compare the results to the the-oretical one. Can you find all the eigenvalues and eigenvectors?

10. Using the fact that the eigendecomposition of a symmetric matrix gives orthogonal matrices (which are easy to invert) propose a method to solve the linear system

x = b

11. Solve the system for b = 0.01^2*sin.(2π*(0.01:0.01:0.99)) and

as a 99 × 99 matrix. Plot the resulting vector x. Plot on the same axis sin.(2π*(0.01:0.01:0.99))/4π^2. Is it similar to x? We have just solved the boundary value problem ″ = sin(2 ) with (0) = (2 ) = 0. In the next problem set we will see how to do this even quicker using Fourier analysis.

Exercise 2: Low-rank approximation

In this problem we will use the SVD to produce low-rank approximations for matrices. We will then use this to compress images and write fast multiplication algorithms.

3

In the last problem set we saw that up algorithms. A matrix is of rank

we could exploit structural zeros to speed if we can write it in the form

= ∑ u v⊤

where the u and v are vectors of length , So u v⊤ is a matrix of size ( × ), of rank 1.

The SVD of a square × matrix exactly writes in this form:

u

u

u

1

−−

v

⊤

−−

u v

v

⊤1

= Σ

⊤

|

|

|

2

2

⊤

v

= [ 1

2 ⋯

]

−−

⊤⋮

−−

= ∑

|

|

|

⋱

=1

−−

−−

Truncating the summation after the largest singular values results in a rank-approximation of the matrix . In fact the Eckhart–Mirsky–Young theorem shows this is the best rank approximation to .

Define the rank approximation of

= ∑ u v⊤

where u and v are the singular vectors of and are the singular values of the matrix.

1. Write a function that takes in a matrix A and uses the SVD to construct its rank- approximation, lowrank_approx(A, r). You may use the svd function from the LinearAlgebra standard library package as U, Σ, V = svd(A).

2. Make a square matrix M of size (20 × 20) that is all zero, except for an axis-aligned cross centred in the centre. Each arm should be of total width 2 and extend a distance 5 out from the centre in each direction parallel to the axes. The cross should consist of 1s.

3. Remembering that the rank is the column rank, i.e. the number of linearly-independent columns, how much should the rank be? Use the SVD of M to confirm this. What does the rank-1 approximation of M look like?

4. Now add small random gaussian noise using the randn function, of inten-sity 0.1. Plot the new matrix (using heatmap). How should this affect the (column) rank of the matrix?

4

5. Plot the singular values as a function of . What do they tell us? What is a suitable rank- approximation to take? What does it look like?

6. Now let’s apply these ideas to a real image. Use the Images.jl package to load a test image using the following code. (Remember that you may need to install the relevant packages):

using Images, TestImages

img = testimage("mandrill")

imgmat = Array(channelview(Float64.(Gray.(img))))

heatmap(imgmat, c=ColorGradient(:grays), showaxis=false, clims=(0,1), yflip=true)

Here imgmat is a standard Julia matrix that contains greyscale pixel val-ues.

7. Plot the rank-50 approximation to imgmat (using heatmap). How does it compare to the original?

8. Plot the singular values as a function of . You should see an “elbow” shape. Approximately where does it occur?

9. Create an interactive visualization that shows the low-rank approximation for different values of . What do you observe? After which are you happy with the quality of the image. Is it related to where the elbow is?

Exercise 3: Dynamic mode decomposition

In this problem we will use the SVD to predict the future!

Suppose that we have some data that satisfies an ODE,

ẋ= f(x)

We can solve the ODE and take time snapshots of the result which we can stack into a matrix as follows:

,

= [

x

|

x

|

x

|

( )

( +1) ⋯

( )]

|

|

|

We will assume that the dynamics is linear. This means that we expect that two consecutive snapshots should satisfy

x +1 = x

for some matrix .

5

If we have snapshots then we can write this as a matrix equation,

2, = 1,( −1)

1. Suppose that we have an eigen-decomposition for , i.e. = Λ −1, and that we can write the initial time snapshot as a sum over the eigenbasis of ,

x0 = ∑ w

where w is the th eigenvector of . Show that future time snapshots can be written as

x = ∑ w

where is the th eigenvalue of .

2. We are now going to try and find the eigendecomposition of without ever constructing it! Start by calculating the SVD of 1,( −1) = Σ ⊤. Find an expression for ⊤ in terms of 2, , Σ, , . We can then calculate the eigenspectrum of A by taking the eigenvalue decomposition of ⊤ since they are similar and using use the result in [1.5].

3. Write a function that calculates the eigenspectrum of given 1, . In-stead of using the full SVD, use a truncated SVD keeping only the first singular values and singular vectors; i.e. use the matrices Vr = V[:, 1:r], Ur = U[:, 1:r] and Σr = Σ[1:r, 1:r] in the expression above for , , Σ. Your function should be of the form aspectrum(X, r). The reason for truncating is in case A is not full rank, in which case some terms in Σ might be nearly 0 and dividing by them would be a bad idea.

4. Test your function on a random 10 × 10 matrix generating some data for of size (10×11) from it starting from a random 0. Compare the re-sults of eigen(A) with aspectrum(X, 10). Remember that eigenvectors are the same upto a multiplicative constant.

5. We are now going to apply this to some dynamical data. Use the ODE integrator code below to generate some data for coupled oscillators with 1(0) = 0.5 and all the others (0) = 0 for N = 10 from = 0 to = 10. The coupled system can be written as

1̈

−2

1

1

1

2̈

1

−2

2

3

=

1

−2 1

3

̈

⋮

⋱ ⋱ ⋱

⋮

̈

1

−2

6

or as a system of first order equations,

1

1̇

0

1

1

1̇

−2

0

1

1

2̇

0

0

0

2

2

1

0

−2

0

1

2

̇

3

3

=

0

0

0

0

1

3̇

1

3

̇

1

0

−2

0

⋮

⋱⋱⋱⋱⋱

⋮

0

0

0

0

1

The output

̇

̇

1

0

−2

0

1, 1̇ , 2,

2̇ , ⋯

is a data matrix with rows

Generate a plot of 1( ) to check that everything went according to plan.

6. Split the data into two parts 1 and 2, the first half from = 0 to = 5 = 5 to = 10. Calculate the spefctrum of withandthesecondhalf

= 10 using 1.

7. Use the first column of 2 as the initial condition. Use the spectrum you found to predict the future dynamics. [Hint: use the initial condition to find the s, which is a matrix solve. Then use the equations in part 1.1 to calculate the prediction.]

8. Plot the prediction for the 10 springs on the same axis as the true solution. What happens?

9. Repeat [3.6–3.7] for = 15 and = 20. What do you observe?

ODE Code:

struct RKMethod{T}

c::Vector{T}

b::Vector{T}

a::Matrix{T}

s::Int

# Make sure that the matrices and vectors have the correct size

function RKMethod(c::Vector{T}, b::Vector{T}, a::Matrix{T}, s::Int) where T lc = length(c); lb = length(b); sa1, sa2 = size(a)

if lc == sa1 == sa2 == s-1 && lb == s

new{T}(c, b, a, s)

else

error("Sizes should be (s = $s) : \n length(c) = s-1 ($lc) \n length(b) = s

end

end

end

7

function (method::RKMethod)(f, x, t, h)

# Extract the parameters

a, b, c, s = method.a, method.b, method.c, method.s

• Vector to hold the k terms k = [f(t, x)]

for i in 1:s-1

tn = t + c[i]*h

xn = x + h*sum(a[i, j] * k[j] for j in 1:i)

push!(k, f(tn, xn))

end

return x + h * sum(b.*k)

end

function integrate(method, f, x0, t0, tf, h)

• Calculate the number of time steps N = ceil(Int, (tf - t0) / h)

hf = tf - (N - 2)*h

#initiate tracking vectors

xs = [copy(x0)]

ts = [t0]

#step

for i in 1:N-1

push!(xs, method(f, xs[i], ts[i], h))

push!(ts, ts[i] + h)

end

# Special last step

push!(xs, method(f, xs[N-1], ts[N-1], hf))

push!(ts, ts[N-1] + hf)

return ts, xs

end

c = [1/2, 1/2, 1]

b = (1/6) .* [1, 2, 2, 1]

a = [1/2 0 0; 0 1/2 0; 0 0 1]

rk4 = RKMethod(c, b, a, 4)

function build_springmat(N)

springmat = zeros(2N,2N)

for i = 1:2N

8

(i < 2N) && (springmat[i, i+1] = 1.0)

if iseven(i)

springmat[i, i-1] = -2.0

(i > 3) && (springmat[i, i-3] = 1.0)

end

end

springmat

end

N=10

const spmat = build_springmat(N)

spf(t, x) = spmat*x

x0 = zeros(2N); x0[1] = 0.5

ts, xs = integrate(rk4, spf, x0, 0.0, 20.0, 0.005); X = hcat(xs...)

plot(ts, X[1:2:2N, :]', leg=:outertopright, box=:on)

Exercise 4: Fourier integrals

In lectures we saw that we could write periodic functions in the form

̂

( ) = ∑

=−

where ̂= 1 ∫2 ( ) − .

2 0

1. Consider the saw-tooth function,

( ) =

mod

0 ≤ <

{2 −

≤ <2

{

( , 2 ))

(

else

Calculate the Fourier series coefficients analytically.

2. Write a function fourier_coefficients(f::Vector, n::Int) that takes in a vector of samples of uniformly distributed over [0, 2 ] and returns an approximation to ̂by calculating the integral using the trapezoidal rule.

3. Now calculate ̂using your trapezoidal code using 100 points and = −3, −2, … , 3. How do they compare to the theoretical results?

4. Fix to be 1. Plot the relative error between the theoretical result and the result from using fourier_coefficients for calculating 1̂using a number of points between 10 and 1000. What does the convergence look like? Does it agree with what we discussed in class?

9

5. Now consider the smooth periodic function exp(cos( )). Repeat [4.3– 4.4] for this function. The analytical result is given by ̂ = | |(1). You can calculate this using the besseli(abs(n), 1) in SpecialFunc-tions.jl.

6. Plot the magnitude of ̂as a function of for the two functions. How do the coefficients decay? Is this what you expected?

10