\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}
  \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 4},
            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 4}

\input{preamble.tex}
\title{Homework}



\begin{document}
\maketitle

Please answer the following questions in complete sentences in submit
the solution on Blackboard October 14th, 2016 by midnight. Note: You can
have an optional extension to October 16th, 2016 by 5pm, but then we
cannot guarantee that you will get your graded homework back before the
midterm.

\textbf{Updates 2016-10-12} * Changed Matlab to Julia in a few places. *
Removed one problem that we won't cover the material (Problem 5.2)

\section{Homework 4}\label{homework-4}

\subsection{Problem 1 (10 points)}\label{problem-1-10-points}

Consider a certain corporation that has three fields of operation; it
mines coal, produces gasoline, and generates electricity. Each of these
activities makes use of varying amounts of three products. Suppose that
in order to produce one unit of coal, the corporation consumes
\[ \bmat{0 \text{ units of coal} \\ 1 \text{ unit of gasoline } \\ 1 \text{ unit of electricity } } \]
to produce one unit of gasoline, the corporation consumes
\[ \bmat{0 \text{ units of coal} \\ 1/5 \text{ unit of gasoline } \\ 2/5 \text{ unit of electricity } } \]
and to produce one unit of electricity,
\[ \bmat{1/5 \text{ units of coal} \\ 2/5 \text{ unit of gasoline } \\ 1/5 \text{ unit of electricity } } \]

Let
\[ \vx = \bmat{x_1 = \text{ number of units of coal produced} \\ x_2 = \text{ number of units of gasoline produced } \\ x_3 = \text{ unumber of units of electricity produced} }. \]
This is called a production vector. Suppose the corporation needs to
produce 100 units of each product to sell above and beyond its internal
requirements.

Derive and solve a linear system for the amount of each material the
corporation has to produce to meet this demand. (Hint: don't forget to
take the corporations own use into account!)

\subsection{Problem 2 (10 points)}\label{problem-2-10-points}

We'll derive a matrix equation for the PageRank linear system. For this
problem, intuitive but \emph{logical} arguments suffice and you need not
formally prove these statements. Although, proofs are the most
\emph{logical} arguments possible and are accepted and encouraged!

Consider the adjacency matrix of a graph \(\mA\). Let \(\mD\) be a
diagonal matrix with the number of out-links from each node on the
diagonal. (Assume that each node has at least one out-link.)

The matrix \(\mD^{-1}\) is also a diagonal matrix where we replace each
diagonal entry with it's reciprocal. (Note, if this isn't a familiar
fact, take a moment to derive it. The property that
\(\mD \mD^{-1} = \mI\) is handy to use.)

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
  (5 points) Argue or prove that if \(\vx\) is a vector giving the
  probability that the surfer is on each page, then
  \(\vy = \mA^T \mD^{-1} \vx\) is a vector giving the probability that
  the surfer is on each page after following one link from the graph.
\item
  (5 points) Now let \(\vx\) be the probability that the surfer is on
  each page after browsing for a really long time. Recall the argument
  from class, and argue or prove. why:
  \[ \vx = \alpha \mA^T \mD^{-1} \vx + (1-\alpha) \ve /n \] where
  \(\ve\) is the vector of all ones and \(n\) is the number of nodes in
  the graph.
\item
  (0 points) Rewrite this into the linear system:
  \[ (\eye - \alpha \mA \mD^{-1}) \vx  = (1-\alpha) \ve /n. \] This give
  us an easy way to compute PageRank scores on a computer without
  simulating the Markov chain!
\end{enumerate}

\subsection{Problem 3 (5 points)}\label{problem-3-5-points}

Consider the LU factorization from the quiz.

\[ \underbrace{\left[ \begin{array}{rrrr}
1 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1 \\
0 & 1 & 0 & 0 \\
\end{array} \right]}_{\mP}
\underbrace{\left[ \begin{array}{rrrr}
1.9410 & -0.3510 & 0.7352 & 0.4583 \\ 
-0.1639 & -1.7494 & 0.0623 & -0.6425 \\ 
0.1269 & 1.8231 & -0.6234 & -2.1405 \\ 
-1.1361 & 0.2914 & -1.4423 & 0.2333 \\ 
\end{array} \right]}_{\mA} =
\] \[
\underbrace{\left[ \begin{array}{rrrr}
1.0000 & 0.0000 & 0.0000 & 0.0000 \\ 
0.0654 & 1.0000 & 0.0000 & 0.0000 \\ 
0.5853 & 0.0465 & 1.0000 & 0.0000 \\ 
-0.0844 & -1.9637 & 0.5331 & 1.0000 \\ 
\end{array} \right]}_{\mL} 
\underbrace{\left[ \begin{array}{rrrr}
1.9410 & -0.3510 & 0.7352 & 0.4583 \\ 
0.0000 & 1.8460 & -0.6715 & -2.1705 \\ 
0.0000 & 0.0000 & -0.9807 & 0.6026 \\ 
0.0000 & 0.0000 & 0.0000 & -3.0168 \\ 
\end{array} \right] }_{\mU}
\]

Report the errors in \(\mL\) and \(\mU\).

\subsection{Problem 4 (The least squares cluster! 25
points)}\label{problem-4-the-least-squares-cluster-25-points}

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
  Chapter 7, problem 14a; but you may not use Julia's ``backslash''
  operation for this part as mentioned in the book. Instead, you must
  write your own least squares solver using the Julia function
  \texttt{qr}. (See page 170 in the text for the algorithm.)
\item
  Chapter 7, problem 14b.
\item
  Chapter 7, problem 14, comment on the differences you see.
\item
  Chapter 7, problem 15b.
\item
  Chapter 7, problem 16; this is the teams problem from class. Solve
  this least-squares problem where Purdue (T1) gets score \(r_1 = 100\).
  And you must derive a least-squares problem where \(r_1\) is not a
  variable. (Hint: it should be a 5-by-3 matrix for the least-squares
  problem.)
\end{enumerate}

\subsection{Problem 5 (Accuracy and conditioning, 10
points)}\label{problem-5-accuracy-and-conditioning-10-points}

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
  Chapter 6, problem 2
\item
  \textbf{Removed, you do not have to complete this problem.} Show than
  an orthogonal matrix has condition number 1. (Hint: the
  submultipliative property is very helpful here.)
\end{enumerate}

\subsection{Problem 6 (Iterative methods, 15
points)}\label{problem-6-iterative-methods-15-points}

This problem asks you to report on and explain the results of the
following Julia script! Your explanations should be simple. Don't think
too much about this and write the first thing that comes to mind.

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
  (1 points) Start by downloading two images.

\begin{verbatim}
download("http://academics.davidson.edu/math/chartier/Numerical/code/dyoung.jpg","dyoung.jpg");
download("http://academics.davidson.edu/math/chartier/Numerical/code/jacobi.jpg","jacobi.jpg");
\end{verbatim}
\item
  (1 points) Show and explain the following output.

\begin{verbatim}
using Images
im1 = load("jacobi.jpg")
im2 = load("dyoung.jpg")

N = max(size(im1)...)

# You don't have to explain what this function is doing
# You just have to describe the result of using it 
# From http://math.mit.edu/~stevenj/18.303/fall13/lecture-10.html
diff1(M) = [ [1.0 zeros(1,M-1)]; diagm(ones(M-1),1) - eye(M) ]
sdiff1(M) = sparse(diff1(M))
# make the discrete -Laplacian in 2d, with Dirichlet boundaries
function Laplacian(Nx, Ny)
   Dx = sdiff1(Nx) 
   Dy = sdiff1(Ny) 
   Ax = Dx' * Dx
   Ay = Dy' * Dy
   return kron(speye(Ny), Ax) + kron(Ay, speye(Nx))
end

# Normalize A to have 1 on the diagonal
# (This simplifies many thing that would 
# otherwise cause us to divide everything
# by 4.)
A = Laplacian(size(im1,1),size(im2,2))
A = A/4

## Show a piece of the matrix A
% A is a sparse matrix, so what does this show?
A[1:10,1:10]  

## Try this one too
full(A[1:10,1:10])
\end{verbatim}
\item
  (1 points) Show and explain the following output.

\begin{verbatim}
# 
using Plots
function myspy(A::SparseMatrixCSC;ms=5.0)
    xn,yn = size(A)
    ai,aj,av = findnz(A)
    scatter(ai,aj,zcolor=av,ms=ms,legend=false,
        yflip=true, aspect_ratio=:equal,
        m=(:auto,1.0,stroke(0)),
        c=:viridis,
        grid=false,
        xlim=(0.5, xn+0.5), ylim=(0.5, yn+0.5),)
end

# You only have to explain this one!
myspy(A)

# And this one
myspy(A;ms=1)
plot!(xlim=(1,512),ylim=(1,512))
\end{verbatim}
\item
  (1 points) Show and explain the following output.

\begin{verbatim}
# What does multiplication by A do?
X = convert(Array{Float64,2},data(im1))'
x0 = reshape(X,N*N)
x = A*x0;
grayim(reshape(4*x,N,N)') # scaling by 4 makes the picture nicer
\end{verbatim}
\item
  (1 points) Show and explain the following output.

\begin{verbatim}
## Create b such that image2 is answer to Ax = b
b = A*reshape(convert(Array{Float64,2},data(im2))', N*N)
x0 = reshape(convert(Array{Float64,2},data(im1))', N*N)

## Show the initial vector
grayim(reshape(x0,N,N)')

## Show the right-hand-side vector
grayim(reshape(4*b,N,N)')
\end{verbatim}
\item
  (5 points) Explain what the linear system \(\mA \vx = \vb\) for this
  problem solves. (Hint: you may want to look at the next two problems
  first, which actually solve it!)
\item
  (1 poins) Run the Richardson method on this problem.

\begin{verbatim}
# Run the Richardson method
# This is incredibly slow in Julia
omega = 1;
x = x0;
niter = 20;

@gif for i=1:niter
    r = b-A*x;
    x = x + omega*r;
    @show i

    heatmap(-reshape(x,N,N);c=:grays,yflip=true,colorbar=false)
end
\end{verbatim}
\item
  (4 points) Modify the above code to answer the following questions.
  (These are simple modifications) What does the solution look like
  after 20 iterations? How many iterations does it take before you think
  the solution ``looks right''? Hint, it may be useful to use the
  following code.

\begin{verbatim}
function richardson(A,b,x0,omega,niter)
    x = x0
    for i=1:niter
        r = b-A*x;
        x = x + omega*r;
    end
    return x
end
omega = 1.0
# This almost works!
x = richardson(A,b,x0,omega,niter)
grayim(reshape(x,N,N)')

# And play it to get the information you need.  
\end{verbatim}
\end{enumerate}

\subsection{\texorpdfstring{Problem \(\star\) (Not graded, but good
practice for the
midterm}{Problem \textbackslash{}star (Not graded, but good practice for the midterm}}\label{problem-star-not-graded-but-good-practice-for-the-midterm}

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
  Chapter 7, problem 6
\item
  Let \(\mA\) be non-singular and let \(\lambda, \vv\) be an eigenvalue,
  eigenvector pair of A. Show that \(1/\lambda, \vv\) are an eigenvalue,
  eigenvector pair of \(\mA^{-1}\).
\item
  ( Show than an orthogonal matrix has condition number 1. (Hint: the
  submultipliative property is very helpful here.)
\end{enumerate}

\end{document}
