a logo for the course

Matrix Computations

David Gleich

Purdue University

Fall 2017

Course number CS-51500

Tuesday and Thursday, 10:30-11:45am

Location Forney B124


Homework 11

Homework 11

Please answer the following questions in complete sentences in a typed, clearly prepared manuscript. This homework is due by Monday, December 4th, 2017, 4:30am. However, everyone is given an automatic, full-points extension until Wednesday, December 6th, 2017 at 4am if you want your homework graded by December 7th in class. This is designed to give you maximum flexibility in order to plan your time. We highly suggest sticking to the original due date, but there is absolutely no penalty for taking more time.

Problem 0: Homework checklist

Problem 1: Using GMRES

It's based on an assignment from Michael Saunders' iterative methods class.

Download the poisson2D.mat file from Tim Davis's sparse matrix repository.
From Julia, you can use the following code to setup and from the linear system.

    using MAT
    file = matread("poisson2D.mat")
    A = file["Problem"]["A"]
    A = (A+A')./2 # make it symmetric
    b = file["Problem"]["b"]
    b = vec(b)

The matrix A is symmetric-positive definite. But we are going to look at solving where . We'll declare a method converged if the relative residual has -norm less than .

  1. Plot the relative residuals for 100 iterations of the MINRES (minres from IterativeSolvers.jl) method and the GMRES (gmres from IterativeSolvers.jl) method restarted every 30 iterations on a semilogy scale. Describe what you observe, in particular, do the methods converge? For convenience, here is some code to help you get started:

    using IterativeSolvers
    x = minres(A,b,tol = mytol, maxiter = mymaxiter)
    x = gmres(A,b,restart = myrestart, tol = mytol, maxiter = mymaxiter)
    
  2. If you use an incomplete cholesky ichol or incomplete LU ilu preconditioner with GMRES, it will converge. How many iterations does it take? Does MINRES benefit from an incomplete Cholesky preconditioner? We provide you with the ichol and ilu code. Download the file hw11_functions.jl from here http://www.cs.purdue.edu/homes/dgleich/cs515-2017/julia/hw11_functions.jl. The IterativeSolvers package still has some rough edges. We will help you set that up. Here is the code to help you get started

    using ItertaiveSolvers
    include("hw11_functions.jl")
    Lchol = ichol(A)
    M = MyPreconditionerType(Lchol)
    x = gmres(A,b,restart = 30, tol = mytol, maxiter = mymaxiter,Pl=M)
    

    For MINRES, it turns out that the minres function in IterativeSolvers still doesn't accept preconditioners. We use another package to test minres.

    # note that the Krylov package is not a registered Julia package 
    # yet but you can still add it as follows
    Pkg.clone("https://github.com/JuliaSmoothOptimizers/Krylov.jl.git")
    using Krylov
    using LinearOperators
    L = ichol(A)
    P1 = opInverse(LowerTriangular(L))
    Mprecond = P1'*P1
    x = Krylov.minres(A,b,M=Mprecond,itmax=mymaxiter)[1]
    
  3. Show the eigenvalues of the matrix before and after preconditioning. Do you observe any clustering?

Problem 2: (Trefethen and Bau, 28.3)

A real-symmetric matrix has an eigenvalue of of multiplicity , while the rest of the eigenvalues are all in absolute value. Describe an algorithm for finding an orthonormal basis of the 8-dimensional eigenspace corresponding to the domininant eigenvalue.