\documentclass[]{article}
\usepackage{lmodern}
\usepackage{amssymb,amsmath}
\usepackage{ifxetex,ifluatex}
\usepackage{fixltx2e} % provides \textsubscript
\ifnum 0\ifxetex 1\fi\ifluatex 1\fi=0 % if pdftex
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\else % if luatex or xelatex
\ifxetex
\usepackage{mathspec}
\makeatletter % undo the wrong changes made by mathspec
\let\RequirePackage\original@RequirePackage
\let\usepackage\RequirePackage
\makeatother
\else
\usepackage{fontspec}
\fi
\defaultfontfeatures{Ligatures=TeX,Scale=MatchLowercase}
\fi
% use upquote if available, for straight quotes in verbatim environments
\IfFileExists{upquote.sty}{\usepackage{upquote}}{}
% use microtype if available
\IfFileExists{microtype.sty}{%
\usepackage{microtype}
\UseMicrotypeSet[protrusion]{basicmath} % disable protrusion for tt fonts
}{}
\usepackage{hyperref}
\hypersetup{unicode=true,
pdftitle={Homework 11},
pdfborder={0 0 0},
breaklinks=true}
\urlstyle{same} % don't use monospace font for urls
\IfFileExists{parskip.sty}{%
\usepackage{parskip}
}{% else
\setlength{\parindent}{0pt}
\setlength{\parskip}{6pt plus 2pt minus 1pt}
}
\setlength{\emergencystretch}{3em} % prevent overfull lines
\providecommand{\tightlist}{%
\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}
\setcounter{secnumdepth}{0}
% Redefines (sub)paragraphs to behave more like sections
\ifx\paragraph\undefined\else
\let\oldparagraph\paragraph
\renewcommand{\paragraph}[1]{\oldparagraph{#1}\mbox{}}
\fi
\ifx\subparagraph\undefined\else
\let\oldsubparagraph\subparagraph
\renewcommand{\subparagraph}[1]{\oldsubparagraph{#1}\mbox{}}
\fi
\title{Homework 11}
\input{preamble.tex}
\title{Homework}
\begin{document}
\maketitle
\section{Homework 11}\label{homework-11}
Please answer the following questions in complete sentences in a clearly
prepared manuscript and submit the solution by the due date on
Blackboard (around Sunday, December 2nd, 2018.) Note that there will be
a full-points extension opportunity.
Remember that this is a graduate class. There may be elements of the
problem statements that require you to fill in appropriate assumptions.
You are also responsible for determining what evidence to include. An
answer alone is rarely sufficient, but neither is an overly verbose
description required. Use your judgement to focus your discussion on the
most interesting pieces. The answer to ``should I include `something' in
my solution?'' will almost always be: Yes, if you think it helps support
your answer.
\subsection{Problem 0: Homework
checklist}\label{problem-0-homework-checklist}
\begin{itemize}
\item
Please identify anyone, whether or not they are in the class, with
whom you discussed your homework. This problem is worth 1 point, but
on a multiplicative scale.
\item
Make sure you have included your source-code and prepared your
solution according to the most recent Piazza note on homework
submissions.
\end{itemize}
\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}
\subsection{Problem 1: Using GMRES with
Preconditioners}\label{problem-1-using-gmres-with-preconditioners}
\emph{For this problem, you must use Julia.}
It's based on an assignment from
\href{http://www.stanford.edu/class/msande318/}{Michael Saunders'
iterative methods class}.
Download
\href{https://www.cs.purdue.edu/homes/dgleich/cs515-2018/homeworks/poisson2D-data.csv}{\texttt{poisson2D-data.csv}}
and
\href{https://www.cs.purdue.edu/homes/dgleich/cs515-2018/homeworks/poisson2D-rhs.csv}{\texttt{poisson2D-rhs.csv}}
and make sure you can load them via the following script. These files
were originally MAT files from \url{https://sparse.tamu.edu/mat/FEMLAB/}
\begin{verbatim}
using DelimitedFiles
using SparseArrays
data = readdlm("poisson2D-data.csv",',')
A = sparse(Int.(data[:,1]), Int.(data[:,2]),(data[:,3]))
A = (A + A')./2
b = vec(readdlm("poisson2D-rhs.csv"))
\end{verbatim}
The matrix A is symmetric-positive definite. But we are going to look at
solving \((\mA - \sigma \mI) \vx = \vb\) where
\(\sigma = 1.7 \times 10^{-2}\). We'll declare a method converged if the
relative residual has \(2\)-norm less than \(10^{-6}\).
\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
Plot the relative residuals for 100 iterations of the MINRES
(\texttt{minres} from \texttt{IterativeSolvers.jl}) method and the
GMRES (\texttt{gmres} from \texttt{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:
\begin{verbatim}
using IterativeSolvers
x, hist_min = minres(A,b,tol = mytol, maxiter = mymaxiter, log=true)
x, hist_gm = gmres(A,b,restart = myrestart, tol = mytol, maxiter = mymaxiter, log=true)
# look at hist_min.data[:resnorm]
\end{verbatim}
\item
If you use an incomplete cholesky \texttt{ichol} or incomplete LU
\texttt{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 \texttt{ichol} and
\texttt{ilu} code. Download the file \texttt{hw11\_functions.jl} from
here
\url{http://www.cs.purdue.edu/homes/dgleich/cs515-2018/homeworks/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
\begin{verbatim}
using IterativeSolvers
include("hw11_functions.jl")
Lchol = ichol(A)
M1 = LowerPrecond(ichol(A))
x, hist_pgm = gmres(A,b,restart = myrestart, tol = mytol, maxiter = mymaxiter,Pl=M, log=true)
\end{verbatim}
For MINRES, it turns out that the \texttt{minres} function in
IterativeSolvers still doesn't accept preconditioners. We use another
package to test \texttt{minres}.
\begin{verbatim}
using Krylov # ] clone https://github.com/JuliaSmoothOptimizers/Krylov.jl.git
using LinearOperators
P1 = opInverse(LowerTriangular(ichol(A)))
Mprecond = P1'*P1;
x, hist_pmin = Krylov.minres(A,b,M=Mprecond,itmax=mymaxiter, rtol=mytol)
\end{verbatim}
\item
Show the eigenvalues of the matrix before and after preconditioning.
Do you observe any clustering?
\end{enumerate}
\subsection{Problem 2: Using Multigrid to solve Poisson's
equation}\label{problem-2-using-multigrid-to-solve-poissons-equation}
\emph{For this problem, if you don't use Julia, you will have to convert
a series of multigrid codes to your own language. This isn't too hard as
it's only about 100 lines of code. However, in one step, you will have
to use a sparse direct solver. This is simple to do in Julia. If you
know how to do it in other languages, then this should be easy. However,
you are on your own!}
Download the codes
\url{http://www.cs.purdue.edu/homes/dgleich/cs515-2018/homeworks/multigrid_functions.jl}
\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
Using these codes, you can solve Poisson's equation as follows:
\begin{verbatim}
X = solve_poisson_direct(poisson_setup(nx,ny, (x,y) -> 1))
\end{verbatim}
and plot a solution via
\begin{verbatim}
using Plots
pyplot()
surface(X)
\end{verbatim}
(which is what a previous homework asked you to do!)
For this problem, we will always have \texttt{n\ =\ nx\ =\ ny}. How
long does it take to solve this for \(n = 31,63,127,255,511,1023\).
(This worked great on my computer in terms of memory usage. But your
milage may vary, if you run out of memory, show the largest problem
you can solve!)
Plot the times on a log-log scale. If you double the problem size,
about how much longer does the computation need?
\item
Explain what the \texttt{apply\_poisson\_jacobi} function does. Where
is the matrix?
\item
Use this method to implement the Jacobi method to solve a linear
system? For a problem with \(n = 31\), estimate how much decrease in
error occurs at each iteration of Jacobi. That is, if \(\vx\itn{k}\)
and \(\vx\itn{k+1}\) are sucessive iterates with error
\(e_k = \normof{\vx\itn{k} - \vx} / \normof{\vx}\) and \(e_{k+1}\),
respectively, then what does the trend in \(e_{k+1}/e_k\) show? (Note,
you probably have to run a few hundred to a few thousand iterations to
see these values start to converge.)
\item
Write a function \texttt{apply\_poisson\_gauss\_seidel} to implement
the Gauss-Seidel without building the matrix. Show the same error
trend for Gauss-Seidel compared with Jacobi.
\item
Estimate the same rate for Jacobi and Gauss Seidel for a problem with
\(n=63\). Given that each iteration of Jacobi and Gauss-Seidel takes
\(O(n^2)\) time, give a rough estimate of the number of iterations
needed to reduce the error to a fixed value of \(\eps\).
\item
Now we get to play around with Multigrid! The function
\texttt{interpolate} takes a \(n = 31\) vector and returns a
\(n = 63\) vector. Use \texttt{interpolate} to evaluate the error that
results from solving a \(n=31\) Poisson problem and interpolating the
solution to a \(n=63\) solution. Show a mesh-plot of the absolute
value of the error and the norm of the error.
\item
What happens when you run a Jacobi iteration after interpolating the
solution? Describe how the error changes in terms of the mesh-plot as
well as the the norm of the error.
\item
Explain what the \texttt{simple\_multigrid} function does. Estimate
the decrease in error for one iteration of the loop for \(n = 31\) and
\(n = 63\).
\item
Explain what the \texttt{apply\_poisson\_multigrid} function does and
how \texttt{poisson\_multigrid\_error} compares with
\texttt{simple\_multigrid}. Estimate the decrease in error for one
iteration of the loop for the same values of \(n\).
\item
In practice, we can't use the error as we won't have the true
solution. Use the \texttt{poisson\_multigrid\_residual} method to
solve Poissons equation with 25 iterations for
\(n = 31,63,127,255,511,1023\). Show the times on a log-log scale
along with the times from the direct method in part 1. Can you solve a
10,000-by-10,000 problem? (This is a big problem, you'll need O(16GB)
of RAM!)
\end{enumerate}
\subsection{Problem 3: Eigenvalue
estimates.}\label{problem-3-eigenvalue-estimates.}
\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
Use the Lanczos method to estimate the top 5 singular values of the
eigenfaces matrix (Homework 7, problem 3). How many iterations do you
need to get a good estimate. Compare this with the power method for
the dominant vector. (Note that you do not need to be concerned with
the vectors for this part.)
\item
Does the loss of orthogonality of the Lanczos vectors impact the
accuracy of the eigenvectors? Describe your investigation here. (A
rough guide, you should present evidence that shows you are familiar
with the ideas in class. This could be 1 paragraph up to around a
page, and could include figures.)
\end{enumerate}
\end{document}